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Abstract 


-  'Direct  numerical  solutions  of  the  three-dimensional 
time-dependent  Navier-Stokes  equations  are  presented  for  the 
evolution  of  three-dimensional  finite-amplitude  disturbances 
of  plane  Poiseuille  and  plane  Couette  flows.  Spectral  methods 
using  Fourier  series  and  Chebyshev  polynomial  series  are  used. 
It  is  found  that  plane  Poiseuille  flow  can  sustain  neutrally 
stable  two-dimensional  finite-amplitude  disturbances  at 
Reynolds  numbers  larger  than  about  2800.  No  neutrally  stable 
two-dimensional  finite  amplitude  disturbances  of  plane 
Couette  flow  were  found. 

Three-dimensional  disturbances  are  shown  to  have  a  strongly 
destabilizing  effect.  It  is  shown  that  finite-amplitude 
disturbances  can  drive  transition  to  turbulence  in  both  plane 
Poiseuille  flow  and  plane  Couette  flow  at  Reynolds  numbers 
of  order  1000.  Details  of  the  resulting  flow  fields  are 
presented.  It  is  also  shown  that  plane  Poiseuille  flow  can  not 
sustain  turbulence  at  Reynolds  numbers  below  about  500.  — ; — 


1.  Introduction 


One  of  the  oldest  unsolved  problems  of  fluid  mechanics 
is  the  theoretical  description  of  the  inception  and  growth 
of  instabilities  in  laminar  shear  flows  that  lead  to  trans¬ 
ition  to  turbulence.  The  behavior  of  small  amplitude  disturb¬ 
ances  on  a  laminar  flow  is  reasonably  well  understood,  but 

understanding  of  the  behavior  of  finite  amplitude  disturbances 
is  in  a  much  less  satisfactory  state.  There  is  as  yet  no 
close  agreement  between  theoretical  and  experimental  studies 
of  transition  flows. 

In  the  laboratory,  Davies  &  White  (1928) ,  Kao  &  Park 
(1970) ,  and  Patel  &  Head  (1969)  have  shown  that  plane  Poiseuille 
flow  is  unstable  to  finite  amplitude  disturbances  at  Reynolds 
numbers  as  low  as  1000  and  that  initially  turbulent  flow 
remains  turbulent  at  slightly  lower  Reynolds  numbers.  Here  the 
Reynolds  number  is  R  =  Uh/v ,  where  U  is  the  maximum 
downstream  velocity,  h  is  the  half-channel  depth  and  v 
is  the  kinematic  viscosity.  On  the  other  hand,  Nishioka, 

Iida  &  Ichikawa  (1975)  performed  experiments  in  a  low- 
turbulence  wind-tunnel  in  which  they  were  able  to  maintain 
laminar  plane  Poiseuille  flow  at  Reynolds  numbers  as  large 
as  8000.  In  order  to  postpone  transition  to  R  =  8000  , 

Nishioka  et  al  had  to  reduce  the  background  turbulence  level 
to  less  than  0.05%.  At  larger  disturbance  levels,  instab¬ 
ilities  were  obtained  at  lower  (subcritical)  Reynolds  numbers. 
The  experiments  of  Nishioka  et  al  were  performed  in  a 
channel  with  aspect  ratio  (ratio  of  width  to  depth  of 


channel)  of  27.4.  At  lower  aspect  ratios, 

the  channel  geometry  may  induce  significant  three-dimen- 
^  sionality  that  may  drive  transition  at  lower  R  .  This  latter 

effect  may  influence  the  results  of  Kao  &  Park  but 
should  not  affect  those  of  Davies  &  White  and  Patel  &  Head. 

I  Thus,  it  appears  that  the  transition  Reynolds  number 

observed  experimentally  depends  on  both  the  spectrum  and 
amplitude  of  the  initial  two  and  three-dimensional  dist¬ 
urbances  to  the  flow.  Typically,  transition  is  observed 
at  Reynolds  numbers  of  about  1000. 

The  experimental  situation  with  regard  to  plane  Couette 
flow  is  far  less  satisfactory.  Reichardt  (1959)  showed  that 
a  turbulent  flow  is  obtained  at  Reynolds  numbers 
(based  on  half-channel  depth  and  wall  velocity)  as  low  as 
750.  Mollo-Christensen  (private  communication)  has  obtained 
similar  results.  These  experimental  findings  may  be  subject  to 
dispute  because  of  end  effects  which  ,  it  seems  ,  are  difficult,  to  remove  . 
In  summary,  the  only  experimental  evidence  available  to 
date  suggests  that  plane  Couette  flow  undergoes  transition 
at  Reynolds  numbers  similar  to  those  of  plane  Poiseuille 
>  flow. 

Some  insight  into  the  mechanism  of  transition  in  planar 
shear  flows  was  given  by  the  pioneering  experiments  of 
Klebanoff,  Tidstrom,  &  Sargent  (1962)  who  studied  the 


* 
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evolution  of  a  controlled  three-dimensional  disturbance 
in  a  laminar  boundary  layer.  They  found  that  prod¬ 
uction  of  longitudinal  vorticity  by  the  three-dimensional 
disturbance  gives  a  secondary  motion  that  creates  local 
inflectional  profiles;  the  resulting  highly  unstable  profiles  lead 
almost  instantaneously  to  turbulent  spots.  The  key  result  obtained  by 
Klebanoff  et  al  is  that  initially  weak  three-dimensional 
disturbances  may  control  the  nonlinear  development  of  the 
flow  and  its  transition  to  turbulence.  In  this  paper  we 
expand  on  this  idea  by  studying  whether  a  similar  effect 
can  control  the  transition  to  turbulence  in  plane  Poiseuille 
and  plane  Couette  flow. 

Let  us  begin  by  reviewing  theoretical  approaches  to 
these  problems.  The  equations  of  motion  are  the  Navier- 
Stokes  equations 

3v(x, t) 

~  -  +  v  (x,  t)  *Vv(x,  t)  =  -  Vp(x,t)  +vV  v(x,t)  (1.1) 

V* v(x, t)  =  0  (1.2) 

where  v(x,t)  =  (u,v,w)  is  the  velocity  field  at  location 
x  -  (x,y,z)  and  time  t,  p(x,t)  is  the  pressure  (divided 
by  density) ,  and  v  is  the  kinematic  viscosity.  To  date  there 
is  no  compelling  evidence  that  the  Navier-Stokes  equations  are 
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in  any  way  inadequate  on  the  space-  and  time-scales  involved 
in  transition  and  turbulence. 

The  flows  discussed  in  the  present  paper  are  confined 
between  rigid  walls  at  z  =  ±1  and  extend  to  infinity  in 
the  horizontal  directions  x,y.  The  boundary  conditions 
at  the  rigid  walls  z  =  ±1  are  that  the  velocity  of  the 
fluid  must  equal  the  velocity  of  the  wall.  In  plane 
Poiseuille  flow,  the  undisturbed  fluid  motion  is  given  by 

v(x,t)  =  (l-z2,0,0)  ,  p(x,t)  =  -2vx  ;  (1.3) 

this  flow  is  driv  n  by  a  pressure  gradient.  In  plane 
Couette  flow,  the  undisturbed  fluid  motion  is  given  by 


v(x,t)  =  (z,0,0)  ,  p  (x,  t)  =  0  ; 


(1.4) 


this  flow  is  driven  by  the  motion  of  the  walls  at  z  =  ±1. 
For  these  flows  the  maximum  velocity  of  the  undisturbed 
flow  is  1,  so  the  Reynolds  number  based  on  half-channel 
width  is 


The  evolution  of  a  small  disturbance  on  a  plane-parallel 
shear  flow  is  governed  by  the  Orr-Sommerfeld  equation,  which  is 
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(1.6) 


(~^y  -a2-g2)2w  =  iR [  (aU-ui)  (-^  -ct2-B2)  w-aU"w] 
dz^  dzz 

with  boundary  conditions 

w  =  w'  =  0  at  z  =  ±1.  (1.7) 

Here  the  unperturbed  velocity  is  (U(z),0,0)  and  the  small 
disturbance  is  assumed  to  have  the  form 

w(x,t)  =  Re[„(2)ei0'x+isy-i“t]  (1.8) 

where  a  and  3  are  the  wave  numbers  in  the  x  and  y 
directions,  respectively,  and  uj  is  the  (complex)  frequency 
of  the  disturbance.  If  a  and  B  are  real  and  lmu>0, 
then  the  small  disturbance  is  linearly  unstable.  On  the 
other  hand  if  all  such  small  disturbances  to  a  plane- 
parallel  shear  flow  have  lmoj<0,  then  the  shear  flow  is 
linearly  stable. 

The  critical  Reynolds  number  R  is  defined  as 

c 

the  lowest  value  of  R  at  which  there  is  any 
solution  of  the  Orr-Sommerfeld  equation  with  Im  to  =  0  . 

For  R  >  R  ,  linearly  unstable  solutions  of  the  Orr- 
c 

Sommerfeld  equation  may  exist.  In  a  unidirectional  plane  - 
parallel  shear  flow  [v(x,y, z) /u(x,y,z)  is  independent  of 
x,y,z  and  w(x,y,z)  =  0  ],  Squire's  theorem  (see  Lin  1955) 
implies  that  if, at  some  Reynolds  number  R  ,  there  exists 
an  unstable  three-dimensional  disturbance  [3/0 

in  (1.8)]  then 
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there  exists  an  unstable  two-dimensional  disturbance- 

18  =0]  at  a  lower  Reynolds  number.  Therefore,  the  mode 

that  becomes  unstable  at  R  must  be  a  two-dimensional 

c 

mode.  [We  emphasize  for  later  reference  that  at  Reynolds 
numbers  larger  than  Rc  ,  the  most  unstable  solution  of 
the  Orr-Sommerfeld  equation  may  be  three-dimensional 
(see  Michael  1961).] 

Asymptotic  analysis  of  the  Orr-Sommerfeld  equation  for 

plane  Poiseuille  flow  using  recently  improved  WKB  techniques 

leads  to  the  estimate  Rc  =5769.7  (Lakin,  Ng  &  Reid  1978); 

earlier  asymptotic  analysis  had  given  Rc~5360  (see  Lin  1955) . 

Direct  numerical  solution  of  the  Orr-Sommerfeld  equation 

gives  Rc~5772.22  (Orszag  1971b).  The  mode  that  becomes 

unstable  at  R  has  wavenumbers  a=l. 02055  and  8=0.  Thus  the 
c 

theory  of  small  amplitude  disturbances  suggests  that  plane 
Poiseuille  flow  is  unstable  only  for  Reynolds  numbers 
greater  than  5772,  in  contrast  with  the  experimental 
observation  of  possible  transition  to  turbulence  at 
Reynolds  numbers  as  low  as  1000. 

In  plane  Couette  flow,  all  numerical  evidence  suggests 
that  all  modes  of  the  Orr-Sommerfeld  equation  are  stable  at 
all  Reynolds  numbers  (Davey  1973) .  The  absence  of  any 
critical  Reynolds  number  Rc  for  plane-Couette  flow  is  in 
conflict  with  the  available  experimental  evidence  that  this 
flow  undergoes  transition  at  modest  Reynolds  numbers. 

Meksyn  &  Stuart  (1951)  suggested  that  finite-amplitude 
nonlinear  effects  may  permit  the  growth  of  disturbances  at 
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subcritical  Reynolds  numbers.  Meksyn  &  Stuart  introduced 
the  so-called  mean  field  equations  in  which  only  the 
interaction  of  the  mean  flow  with  the  primary  disturbance 
wave  is  retained  and  higher  harmonics  are  neglected.  They 
found  that  finite-amplitude  two-dimensional  disturbances 
to  plane  Poiseuille  flow  are  unstable  at  Reynolds  numbers 
larger  than  about  2900  with  a  threshold  amplitude  of  about 
8%  of  the  centerline  velocity.  Numerical  work  by  Grohne 
(1969)  solving  the  mean  field  equations  gave  a  critical 
Reynolds  number  of  about  2500  (based  on  the  perturbed 
centerline  velocity) . 

Stuart  (1960)  and  Watson  (1960,  1962)  extended  the 

Meksyn-Stuart  theory  to  include  the  second  harmonic  of  the 

two-dimensional  primary  wave.  Their  results  for  plane 

Poiseuille  flow  are  close  to  those  of  Meksyn  &  Stuart. 

Further  work  using  the  Stuart-Watson  method  by  Reynolds  & 

Ratter  (1967)  has  conf  irnedthat  velocity  fluctuations  of  a 

few  percent  can  drive  two-dimensional  finite-amplitude 

instabilities  at  Reynolds  numbers  above  2900.  The  original 

Stuart-Watson  method  involves  expansions  in  both  small  amplitude 

and  small  Im(u>)  about  the  neutrally  stable  modes  at  R  . 

Recently,  Itoh  (1977)  has  extended  the  method  of  Eckhaus 

(1965)  in  order  to  reformulate  the  theory  in  a  manner 

that  he  claims  avoids  the  restriction  to  the  neighborhood  of 

R  ;  a  critique  of  Itoh's  theory  is  given  by  Davey  (1978) . 
c 
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A  slightly  different  approach  to  the  two-dimensional  finite- 
amplitude  instability  problem  has  been  given  by  Pekeris  & 

Shkoller  (1967,1969,1971),  Zahn,  Toomre,  Spiegel  &  Gough  (1974), 
and  Herbert  (1976,1977).  In  these  investigations,  the  two- 
dimensional  Navier-Stokes  equations  are  solved  using  periodic 
boundary  conditions  in  x  by  expanding  the  solution  in  a 
highly  truncated  Fourier  series  in  x.  Pekeris  &  Shkoller 
(1971)  found  two  dimensional  instabilities  of  plane  Poiseuille 
flow  at  Reynolds  numbers  as  low  as  1000;  at  R  =  1000,  instab¬ 
ility  is  achieved  with  a  6%  perturbation,  while  at  R  =  3000 
instability  is  achieved  with  a  0.4%  perturbation.  art 

(1971)  faults  the  accuracy  of  the  Pekeris-Shkol?  (1967,1969) 
calculations.  Our  numerical  solutions  of  the  Na\ '  _-Stokes 
equations  reported  in  Sec.  3  give  no  evidence  of  two-dimensional 
finite-amplitude  instabilities  in  the  vicinity  of  those  predicted 
by  Pekeris  &  Shkoller  (1971)  on  the  basis  of  nonlinear 
stability  calculations  using  expansions  in  eigenfunctions  of 
the  Orr-Sommerfeld  equation. 

Zahn  et  al  solved  the  Navier-Stokes  equations  for  plane 
Poiseuille  flow  by  retaining  only  two  Fourier  modes  in  x 
and  using  an  unequally  spaced  finite  difference  grid  in  z. 

They  found  a  minimum  critical  Reynolds  number  for  two-dimensional 
finite- amplitude  instability  of  2707.  The  instability  at  this 
Reyholds  number  is  achieved  with  a  value  of  a  =  1.3126. 

Herbert  (1976,1977)  performed  a  similar  calculation  using 
up  to  eight  Fourier  harmonics  in  x  and  41  Chebyshev  polynomials 
in  the  z  direction.  He  found  a  critical  Reynolds  number  of 
2935  with  a  corresponding  a  =  1.3231.  Herbert's  calculations 


are  directly  comparable  to  ours  in  that  he  used  essentially 
the  same  numerical  technique  as  we  do.  The  principal  diff¬ 
erences  are  that;  (i)  Herbert  seeks  neutrally  stable  finite- 
amplitude  two-dimensional  modes  by  solving  time-independent 
equations  while  we  solve  the  time-dependent  problemjand 
(ii)  Herbert  uses  up  to  eight  Fourier  harmonics  in  x 
while  we  use  up  to  32.  Our  two-dimensional  calculations  are 
in  good  agreement  with  those  of  Herbert. 

In  summary,  the  best  available  evidence  to  date  suggests 
that  two-dimensional  disturbances  are  unstable  only  for  Reynolds 
numbers  larger  than  about  2800.  Our  numerical  solutions  reported 
in  Sec.  3  confirm  this  result. 

Direct  numerical  calculations  of  the  two-dimensional 

Navier-Stokes  equation  were  performed  by  George,  Heliums 
&  Martin  (1974) .  They  obtained  instability  only  for 
Reynolds  numbers  larger  than  about  3500.  Further  discussion 
of  their  results  is  given  in  Sec.  3. 

The  current  state  of  understanding  of  the  effect  of 
three-dimensional  disturbances  on  plane  Poiseuille  flow 
is  less  well  settled.  Meksyn  (1964)  applied  the  mean 
field  equations  for  plane  Pouiselle  flow  and  found  three- 
dimensional  finite-amplitude  instability  at  R  =  1260, 
but  he  also  found  two-dimensional  instability  at  R  =  1270. 

The  inconsistency  of  these  results  with  other  numerical 
calculations  of  the  mean  field  equations  and  the  Stuart- 
Watson  equations  remains  unexplained,  but,  in  any  case, 
they  do  not  show  a  large  effect  of  three-dimensionality. 
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Zahn  et  al  (1974)  examined  some  three-dimensional  modes 

and  found  them  to  be  at  least  as  stable  as  the  two- 

dimensional  ones.  Most  recently,  Itoh  (1978)  has  extended 

the  Stuart-Watson-Eckhaus  theory  to  three-dimensional  disturbances 

for  Reynolds  numbers  close  to  the  linear  stability  limit 

of  5772.  itoh  found  that  three-dimensional  disturbances 

are  strongly  destabilizing,  but  Davey 1 s  (1978)  criticism  may 

still  apply;  in  any  case,  Itoh's  theory  does  not  apply 
to  strongly  subcritical  Reynolds  numbers.  Also, 

Hocking,  Stewartson  &  Stuart  (1972)  and  Davey, 

Hocking  &  Stewartson  (1974)  studied  the  evolution  of  three- 

dimensional  finite  amplitude  disturbances  to  plane  Pouiselle 
flow  at  supercritical  Reynolds  numbers.  They  found  that  three- 
dimensional  effects  are  destabilizing  above  the  linear  stability 
limit. 

Theoretical  investigations  of  the  finite-amplitude 
stability  of  plane  Couette  flow  are  also  incomplete  at  this 
time.  Kuwabara  (1967)  applied  the  mean  field  equations  of 
Meksyn  &  Stuart  and  found  that  the  minimum  critical  Reynolds 
number  for  finite-amplitude  instability  of  plane  Couette  flow 
to  two-dimensional  disturbance  to  be  R  =  45212  with  the 
unstable  mode  having  wavevector  a  =  13.565,  0=0.  Ellingsen, 

Gjevik  &  Palm  (1970)  Davey  &  Nguyen  (1971),  and  Coffee  (1977) 
used  the  Stuart-Watson  method  to  study  two-dimensional  finite- 
amplitude  instability.  They  found  instability  to  low-amplitude 
disturbances  down  to  Reynolds  numbers  of  1000  and  below.  However, 
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the  absence  of  neutrally  stable  linear  eigenmodes  of  plane 
Couette  flow  casts  some  doubt  on  the  applicability  of  the 
Stuart-Watson  method  which  involves  an  expansion  about  neutral 
stability  (Rcsenblat  &  Davis  1978) .  Lessen  &  Cheifetz  (1975) 
also  studied  the  nonlinear  evolution  of  two-dimensional 
disturbances  of  plane  Couette  flow.  Their  calculations  cast 
doubt  that  any  unstable  two-dimensional  disturbances  exist. 

Herbert  (1977)  reports  inability  to  find  neutrally  stable 
finite  amplitude  solutions  of  plane  Couette  flow  using  highly 
truncated  Fourier  expansions  in  x.  Our  numerical  solutions 
of  the  Navier-Stokes  equations  reported  in  Sec.  3  give  no 
evidence  yet  of  two-dimensional  finite-amplitude  instabilities 
in  the  neighborhood  of  those  predicted  by  Ellingsen  et  al , 

Davey  &  Nguyen,  and  Coffee.  In  summary,  there  are  some 
significant  disagreements  on  the  existence  and  strength  of 
two-dimensional  finite  amplitude  instabilities  of  plane 
Couette  flow. 

In  this  paper,  we  solve  the  Navier-Stokes  equations 
numerically  to  study  quantitatively  the  instability  and 
transition  to  turbulence  of  plane  Poiseuille  and  plane 
Couette  flows.  Since  the  resulting  turbulence  is  strongly 
three-dimensional  and  since  two-dimensional  nonlinear 
disturbances  of  these  laminar  shear  flows  do  not  seem  to 
be  able  to  explain  observed  experimental  results,  we  con¬ 
centrate  on  the  study  of  possible  three-dimensional  mechanisms. 

Some  insight  is  given  by  results  for  transition  in  boundary 
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layer  flows  where  the  experiments  of  Klebanoff  et  al  (1962) 
and  the  theory  of  Benney  and  Lin  (Benney  &  Lin  1960,  Benney 
1961,  1964)  suggest  that  the  secondary  motions  produced  by 
the  interaction  of  three-dimensional  modes  with  two-dimensional 
modes  can  produce  velocity  profiles  that  are  highly  inflectional 
and  unstable.  Numerical  calculations  of  the  Benney-Lin 
equations  by  Antar  &  Collins  (1975)  have  shown  that  this  kind  of 
theoretical  approach  is  in  some  agreement  with  experiments. 
Direct  numerical  calculations  of  the  three-dimensional 
Navier-Stokes  equation  for  boundary-layer  flow  have  shown 
quantitatively  the  strength  of  these  three-dimensional 
effects  in  producing  transition  (Orszag  1976) .  On  the  other 
hand,  direct  numerical  calculations  of  the  two-dimensional 
Navier-Stokes  equations  for  boundary  layer  flow  (Fasel  1976, 
Fasel,  Bestek  &  Schefenacker  1977,  Murdock  1977,  Murdock  & 

Taylor  1977)  do  not  exhibit  explosively  strong  physical 
instabilities  and  the  small  scale  excitation  and  apparent 
randomness  characteristic  of  transition. 

Because  of  the  limited  spatial  resolution  of  our  calculations, 
we  do  not  address  in  detail  the  nature  of  the  small  scale 
flow  structures  that  result  from  flow  breakdown.  Our  goal 
is  to  explain  the  mechanisms  by  which  flows  develop  that  break 
down  to  turbulence.  The  development  of  very  small  scale 
structures  in  these  flows,  as  studied  by  Landahl 
(1972),  is  beyond  the  scope  of  the  present  work. 

In  Sec.  2,  we  discuss  briefly  the  numerical  methods 
used  in  the  present  study.  Then,  in  Sec.  3,  we  present 


results  for  two-dimensional  linear  and  nonlinear  disturbances 
of  plane  Poiseuille  and  plane  Couette  flow.  In  Sec.  4,  we 
present  results  of  calculations  of  three-dimensional  finite- 
amplitude  instabilities  of  these  flews  and  the  resulting 
transition  to  turbulence.  Finally,  in  Sec.  5,  we  summarize 
our  results. 
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2.  Numerical  Methods 

We  solve  the  Navier-Stokes  equations  (1.1) -(1.2) 
expressed  in  rotation  form 

3v(X/t)  2 

- -  =  v(x,  t)  (x,  t)  -  VII(x,t)  +  vV  v(x,t)  ,  (2.1) 

where  w(x,t)  =  V*v(x,t)  is  the  vorticity  and 

1  2 

n  (x,t)  =  p(x,t)  +  2-|v(x,t)|  is  the  pressure  head. 

The  flow  is  assumed  to  take  place  in  the  three-dimensional 
box  0£x£X,  -  Y  y  £  ,  -  1  <  z  <  1 .  At 

z  =  ±1  we  impose  the  boundary  conditions  that  the  fluid 
velocity  match  the  wall  velocity.  In  the  horizontal 
directions,  we  impose  periodic  boundary  conditions  so  that 

v(x+mX,  y+nY,  z,t)  =  v(x,y,z,t)  (2.2) 

for  all  integers  m,  n.  These  periodic  boundary  conditions 
are  consistent  with  the  Navier-Stokes  equations  and  the 
laminar  solutions  (1.3)  for  Poiseuille  flow  and  (1.4)  for 
Couette  flow. 

The  choice  of  periodic  boundary  conditions  in  horizontal 
planes  does  cause  some  problem  with  respect  to  comparison 
with  experiment  since  these  boundary  conditions  are  not 
realized  in  the  laboratory.  There  are  two  justifications 
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for  their  use:  (a)  The  instabilities  of  laminar  flows 
th&t  lead  to  turbulence  are  of  small  spatial  scale  so  that 
boundary  conditions  should  have  little  effect;  and  (b) 
the  spatial  growth  of  a  disturbance  in  a  laboratory 
coordinate  frame  appears  in  an  advected  coordinate  frame 
as  temporal  growth,  similar  to  that  observed  with  the 
boundary  conditions  (2.2).  Transition  experiments  in  a 
flat  plate  boundary  layer  have  been  performed  (Orszag  1976) 
with  proper  inflow-outflow  boundary  conditions  applied  and 
the  results  are  qualitatively  the  same  as  in  the  present 
work.  Also,  Fasel  et  al  (1977)  have  used  inflow-outflow 
boundary  conditions  in  their  numerical  simulations  of  two- 
dimensional  disturbances  to  plane  Poiseuille  flow  with  results 
similar  to  those  obtained  using  periodic  boundary  conditions. 
We  offer  no  more  excuses  for  the  boundary  conditions  (2.2) 
used  here  and  urge  further  study  of  their  effect  by  future 
investigators . 


We  solve  the  Navier-Stokes  equations  in  Eulerian  coor¬ 
dinates  using  the  pseudospectral  method  suggested  by 
Orszag  (1971a) .  An  introduction  to  the  theory  of  spectral 
methods  is  given  in  the  monograph  by  Gottlieb  &  Orszag 
(1977) .  Here  we  summarize  the  implementation  of  spectral 
methods  for  the  present  channel  flow  simulations. 

The  velocity  field  is  represented  using  Fourier  series 
in  x  and  y  and  a  Chebyshev  polynomial  series  expansion 
in  z  .  Thus,  the  velocity  field  is  represented  as 


v(x,t)  = 


|m|  <M 


A 

I  n!  <N 


A  u(m,n,p,t)e2lTi(mx/X+ny/Y)T  (z) 

p=0  '■  p 
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(2.3) 


where  m,  n,  p  are  integers  and  T^(z)  =  cos (p  cos  ^  z) 
is  the  Chebyshev  polynomial  of  degree  p  . 

Equations  for  the  spectral  components  u(m,n,p,t)  are 
obtained  using  a  pseudospectral  method  to  evaluate  the  non¬ 
linear  terms  of  the  Navier-Stokes  equation.  Thus,  the 
rotation  term  v  x  oi  in  (2.1)  is  evaluated  as 


v  x  U  (Xj  ,yR,z  =  v(x..  yk,z  &)  *  oj  (x..  ,yR,  2  J 
(0  <  j  <  2M,  0  <  k  <  2N,  0  <:  S,  <  P 


(2.4) 


where  the  collocation  points  x  ^ ,  y^,  z^  are 

=  jX/(2M),  yk  »  (k-N)Y/(2N),  z  %  *  cosM/P  . 


The  values  of  v(x.,y.,z  )  are  obtained  from  (2.3)  using 
*v  3  ^  ^ 

fast  Fourier  transform  algorithms  [improved  to  take 
advantage  of  the  reality  of  v  and  the  cosine  (Chebyshev) 
transform  in  z  as  described  in  the  appendix  to  Orszag 
(1971a)].  Similarly,  u(x.,y.,z  )  is  evaluated  using 
fast  Fourier  transforms  applied  to  the  curl  of  (2.3); 
for  this  purpose  it  is  helpful  to  note  that 


3y 
3  z 


T 

|m|  <M 


J]  u(1Nm,n,p)e2’i(mx/X+ny/Y,TpU> 

|n|<N  p=0 

(2.5) 
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where 


c 


u 


(1) 


p-1 


u 


(1) 


is  given  in  terms  of  u  by 
(m,n,p-l)  -  u(1>  (m,n,p+l)  = 


the  recurrence  relation 


2p  u(m,n,p)  (l£p£P) ; 


(2.6) 


where  c  =  1  if  q>l  and  cn  =  2,  and  (m,n,P+l)= 

q  —  0  - 

(m,n,P)  =  0  for  all  m,n  .  Also,  we  apply  a  special 
circular  truncation  to  the  spectral  representation  of  v  *  oi 
in  the  x-y  plane  in  order  to  minimize  aliasing  effects 
(Orszag  1971a,  Sec. 6) . 

The  evaluation  of  v  x  us  by  this  algorithm  requires 
9  Fourier  transforms  on  2 MNP  complex  data  points  [three 
transforms  each  to  get  v  in  physical  space,  to  in  physical 
space,  and  v  x  to  back  in  transform  (Fourier-Chebyshev) 
space  (which  is  the  resident  representation  through  most  of 
our  computer  code)]  .  With  M=N=16  and  P  =  32,  each  component 
of  the  velocity  field  is  represented  by  33,792  real  degrees 
of  freedom  (before  the  circular  truncation  in  the  x-y  plane 
is  applied)  ;  evaluation  of  v  x  to  by  the  above  pseu- 
dospectral  algorithm  requires  about  2.5s  on  the  CDC  7600.  In 
contrast,  direct  evaluation  of  the  convolution-like  sums 
that  would  be  obtained  by  formulating  equations  for 
3u(m,n,p,t)/9t  using  a  Galerkin  approximation  procedure 
would  require  about  1000  times  more  computer  time.  It  is 
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noteworthy  that  of  this  speedup  by  a  factor  1000,  the  fast 
Fourier  transform  contributes  only  roughly  a  factor  2;  most 
of  the  speedup  is  due  to  the  reorganization  of  the  calculation 
in  terms  of  transforms  (which  factor  into  a  sequence  of  one¬ 
dimensional  transforms)  be  they  fast  or  not.  This  latter 
result  and  the  result  noted  by  Gottlieb  &  Orszag  (1977) 
that  Legendre  polynomial  expansions  are  especially  accurate 
may  suggest  that  Legendre  polynomials  be  used  in  place  of 
Chebyshev  polynomials  in  future  studies.  We  also  note  that 
one  of  us  (S.A.O)  has  recently  developed  a  'fast*  Legendre 
transform  that  requires  0  (N  Uog2  N)  2/  £og  2  £og2N)  operations 
to  transform  an  N  term  Legendre  polynomial  or  surface 
harmonic  series  to  or  from  physical  space. 

Four  additional  features  of  our  numerical  methods 
deserve  comment:  time  stepping,  determination  of  the 
pressure,  viscous  dissipation,  and  initial  conditions.  A 
fractional  step  method  is  used  in  time  so  that  we  shall 
first  discuss  a  time  step  of  the  inviscid  equation,  then 
the  imposition  of  the  incompressibility  constraint, 
then  viscous  effects,  and  finally  initial  conditions. 

In  order  that  time  steps  not  be  unduly  restricted 
by  convective  stability  conditions  due  to  the  relatively 
la^ge  unperturbed  (laminar)  flow  we  use  a  semi-implicit 
scheme  in  which  the  largest  effects  of  this  parallel  flow 
are  handled  implicitly.  If  the  (x-component  of  the) 
unperturbed  parallel  flow  is  denoted  by  U(z),  then  the  terms 
in  (1.1)  responsible  for  convective  stability  restrictions 
that  trace  to  U(z)  are  just  U(z)8v/3x.  This  term  is  best 


evaluated  in  the  mixed  spectral-physical  space  representation 
in  which  x  and  y  are  represented  as  Fourier  modes  while 
z  is  kept  in  physical  space;  in  this  representation 
U(z)3v/3x  i  ( 2frm/x)  U (z)  u (m,n,  z) /  which  is  diagonalized. 

With  the  semi-implicit  treatment  of  U(z)  ,  it  is  a 

little  tricky  to  obtain  an  efficient  time-stepping  procedure 

t 

for  (2.1)  and  (1.2)  that  is  high-order  accurate  in  time. 

We  use  an  Adams-Bashforth-Crank-Nicholson  (ABCN)  method 

2 

with  global  error  of  order  0(At  )  +  0 (vAt) ,  where  we 
tolerate  first-order  accuracy  on  the  viscous  terms  because 
our  calculations  are  done  with  very  small  values  of  v.  if 
v  is  large,  the  semi-implicit  treatment  of  U(z)  should  be 
sacrificed  for  better  accuracy  in  time.  The  ABCN  scheme 
for  the  nonlinear  terms  in  (2.1)  is  derived  by  writing 
3v/3t  +  U  3v/3x  =  v  x  a)  +  U  3v/3x  and  applying  the  Crank- 
Nicolson  scheme  on  che  left  and  the  second-order  A dams- 
Bashforth  scheme  on  the  right;  the  result  is 


j'.n+l  n 
u  -u 

At 


1  ^rm 

2  X 


U(z) [Qn+1-2fin+Qn_1] 


(2.7' 


in  the  mixed  spectral-physical  space  representation.  Here 
F  =  v  x  u>+  f  in  the  mixed  representation  where  we  have 
generalized  the  problem  slightly  by  introducing  an  external 
force  (mean  pressure  gradient)  denoted  by  f;  the  superscripts 
label  the  time  step  and  no  boundary  conditions  (except 


t  Orszag  &  Deville  (to  be  published)  have  recently  shown  how 
to  obtain  unconditionally  stable,  easily  implementable 
schemes  of  arbitrary  order  accuracy  in  time  for  the  Navier- 
Stokes  equations.  Comparisons  with  the  results  reported 

Present  paper  show  no  appreciable  errors  due  to  the 
lower  order  scheme  used  here. 
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t 


> 


♦ 


» 


t 


♦ 


t 


♦ 


t 


♦ 


periodicity  are  applied  yet.  Eq.  (2.7)  is  formulated 
for  constant  time  steps  At;  for  nonconstant  time  steps 
[as  for  the  first  few  time  steps  (we  use  a  slow  start 
to  minimize  initial  truncation  errors) ] ,  the  coefficients  in 

(2.7)  are  different.  It  is  important  to  emphasize  that 

the  terms  multiplying  U(z)  on  the  right  side  of  (2.7) 

are  all  intermediate  results  at  time  steps  n-1,  n,  and 

n+l;  if  un~  and  un  are  replaced  by  un_  and  u 
3  n  1 

or  ~  un  +  iu  (as  would  be  the  case  if  naive  application 
of  the  Adams-Bashforth  method  were  made)  then  the  scheme 
(2.7)  would  lead  to  errors  of  order  0(At)  after 
application  of  the  pressure  terms  in  (2.1) . 

Once  the  fractional  step  (2.7)  is  made,  it  is 
necessary  to  include  the  effect  of  pressure.  This  is 
done  by  the  fractional  step 


"n+l  "n+l 
v  -  v 


vn 


n+l 


(2.8) 


V«  vn+1  =  0 


(2.9) 


where  vn+1  has  global  error  of  order  0(At2)  despite  the 
fact  that  AtH  has  global  error  0(At)  relative  to  the 

pressure  head  II.  Eqs.  (2. 8) -(2. 9)  are  solved  with  the 
boundary  conditions 

A 

wn+1  =  0,  z  =  ±1,  (2.10) 


21 


> 


» 


I 


♦ 


♦ 
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assuming  no  normal  wall  motion. 

Eqs.  (2.8)-(2,10)  are  best  solved  in  the  full 
spectral  representation  using  a  spectral  tau  method  in 
z  (Gottlieb  &  Orszag  1977,  Sects.  10,  13).  The  result- 

A 

ing  equations  for  u  (dropping  the  time-step  label)  are, 
in  the  full  spectral  representation  given  by 


u(m,n,p)  =  u(m,n,p)  -  iamn(m,n,p)  (0<p<P)  (2.11) 

A  A  A 

v(m,n,p)  =  v (m, n,p)  -  iBnII(m,n,p)  (0<p<P)  (2.12) 


ft 

w(m,n,p)  *  w(m,n,p)  - 

n(1)  (m,n,p)  (0<p<P-2) 

'V  A 

A  A 

iamu(m,n,p)  +  iBnv(m,n,p)  + 

w^(m,n,p)=0  ( 0 <p<P ) 

P  P 

A.  _  A 

l  w(m,n,p)  =  l  (-l)pw(m,n,p)  =  0 

p=0  p=0 


(2.13) 

(2.14) 


(2.15) 


for  all  retained  m, n  ( \ m j <M, in ] <N) ,  where + 

A 

a=  2tt  /X  ,  B-  2v  /Y  and  H  ^  and  w^  are  related 

A  A 

to  H  and  w  as  in  (2.6). 

2  2 

*  If  m  +n  7*  0,  an  efficient  procedure  to  solve 
(2.11)- (2.15)  is  to  rewrite  ( 2 . 11) - (2 . 15)  as  the  nearly 
tridiagonal  system  (see  Sec.  10  of  Gottlieb  &  Orszag  1977) 


t  Note  that  these  definitions  of  a  ,  B  are  consistent  with  a  •  & 
used  in  (1.8)  and  ( 2 . 23) - ( 2 . 27)  below  if  the  linear 
mode  in  (1.8)  is  the  fundamental  mode  in  the  box 
°<x<X,-  4-  Y < y <  j-Y . 
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4p^rr  w (m, n/P-2)  -  (1  +  -L|±l-)w(,n,p)  +  (m, n,P+2) 

^ IP  “A/ 


=  F(m,n,p)  (2^p<P) 


(2.16) 


where  c  =  2,  c  =  1  (p>_l)  ,  e  =  l(p<P),  e  =  0(p>P) 

^  P*  r'  *r 


2  2  2  2 

Y  =  am  +  p  n  ,  and 


c  0f 


e_  ,  „f  e  ..f 


+  <2^p)  (2-17> 


with 


/s  (1)  ^  (1)  ~ 

f  =  -  iamu  (m,n,p)  -  ignv  (m,n,p)  -Yw(m,n,p) 

r 


(2.18) 


A (i)  A (i)  ~ 

and  u  and  v  related  to  u  and  v  by  (2.6) . 

A 

The  system  ( 2 . 15) - (2 . 16)  is  solved  for  w  by  standard 
techniques  in  roughly  the  same  number  of  operations  required 

/v 

to  solve  pentadiagonal  systems.  The  equations  for  w  in 
this  form  are  essentially  diagonally  dominant  so  no  appreciable 
accumulation  of ’roundoff  errors  occurs. 

A  A 

Once  w  is  found  then  II  is  found  from 


A  ~  ~  ft  (1) 

Yll(m,n,p)  =  -  iau(m,n,p) -i8nv(m,n/P)  -  w  (m,n,p) 


(0<p<P) 


(2.19) 
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where  a,{$,Y  are  given  as  before  in  terms  of  m,n,  Then 

A.  A 

u  and  v  are  found  directly  from  (2 .11) - (2 .12) ,  completing 

A 

the  solution  for  u  . 

If  m  =  n  =  0  ,  then  ( 2 . 11) - (2 . 15)  are  easily  solved 
by  appropriate  applications  of  recurrences  like  (2.6). 

Finally,  we  use  a  full  implicit  fractional  step 
to  impose  the  viscous  terms  and  the  boundary  conditions: 


yn+l_~n+l  =  vV2vn+l 


(2.20) 


vn+1  =  V  (z  =  ±1)  (2.21) 

~  ± 

where  V+  is  the  wall  velocity  at  z  =  ±1,  respectively. 

The  full  implicit  scheme  ( 2 . 20) - ( 2 . 21)  is  unconditionally 
stable  but  it  does  induce  global  errors  of  order  0(v^t) . 

Eqs.  (2. 20) -  (2. 21)  are  solved  efficiently  in  a  full  spectral 
representation  using  a  spectral-tau  method;  the  resulting 
equations  are  essentially  tridiagonal  in  the  Chebyshev  index 
and  diagonal  in  the  Fourier  indices  (see  Sec.  10  of  Gottlieb 
&  Orszag  1977)  .  This  completes  the  description  of  the  main 
part  of  our  computer  code. 

Initial  conditions  for  the  runs  reported  in  Sects. 

3-4  usually  consist  of  the  unperturbed  flow  (1.3)  or  (1.4) 
on  which  we  superpose  finite-amplitude  combinations  of 
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two-  and  three-dimensional  eigenmodes  of  the  Orr-Sommerfeld 
equation  (1.6) -(1.8).  The  Orr-Sommerfeld  equation  is 
solved  by  expanding  w(z)  in  the  Chebyshev  series 
(Orszag  1971b) 

P 

w(z)  =  l  a  T  (z) ,  (2.22) 

p=0  p  p 

constructing  equations  for  the  expansion  coefficients 
Sp  by  matrix  methods  (Metcalfe  1974)  ,  finding  the  eigen¬ 
value  by  either  a  global  eigenvalue  routine  based  on 
a  matrix  QR  eigenvalue  analysis  or  a  local  eigenvalue 
routine  based  on  an  inverse  iteration-Rayleigh  quotient 
method,  and  finally  obtaining  the  eigenfunction  by  an 
inverse  iteration  method.  The  accuracy  of  our  eigenvalues 
and  eigenfunctions  is  better  than  1  part  in  10**. 

Once  the  complex  z-velocity  w(z)  of  an  eigenmode 
of  the  Orr-Sommerfeld  equation  is  obtained,  the  x-  and 
y-  velocity  components,  u(z)  and  v(z)  ,  respectively, 
are  obtained  as  follows.  The  z-component  of  the  perturbed  vorticity, 
given  by  n (z) exp [iax+igy-iwt]  with 

T]  (z)  =  iPu(/>  -  iav(z)  ,  (2.23) 
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satisfies  the  linear  inhomogeneous  equation 


a  2  2 

( — j  ~  a  ~  e  )H(z)=  iR[  (aU(z)-u)n  (z)+i8U*  (z)w(z)  ]  (2.24) 

dz 


n  (±1)  “  0  /  (2.25) 

where  R  =  1/v.  Using  the  incompressibility  condition 
(1.2),  it  follows  that 

(a2+82)u(z)  =  -i8n  (z)  +iaw'(z)  (2.26) 

(a2+82)v(z)  =  ian(z)  + iBw* (z) .  (2.27) 

For  two-dimensional  disturbances,  3=0  and  n(z)  =  v(z)  =0. 

In  concluding  this  section  we  summarize  some  features 
of  our  computer  code.  Fourier  series  representations  are 
used  in  x  and  y  and  Chebyshev  series  representations 
are  used  in  z  .  Pseudospectral  methods  are  applied  to  the 
nonlinear  terms  while  tau  methods  are  ipplied  to  the 
pressure  and  viscous  terms.  The  resulting  scheme  is  infinite- 
order  accurate  in  space,  has  no  phase  errors  in  x  and  y  , 
can  resolve  boundary  layers  of  thickness  0(1/P  )  in  z  , 

[the  collocation  points  z^  and  z^  ^  in  (2.4)  are 
located  a  distance  tt  /2P  from  the  walls  at  z  =  ±1] 
and  accurately  imposes  the  boundary  conditions  at  the  rigid 
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walls  and  the  incompressibility  condition  throughout 

the  layer  (see  Gottlieb  &  Orszag  1977  for  discussion  of 

these  features) .  A  fractional  time  step  method  is  used 

2 

with  global  error  of  order  0 (At  +vAt) ;  time  steps  are 
formally  restricted  only  by  convective  stability  rest¬ 
rictions  due  to  the  perturbed  velocity  [although  in 
practice  we  do  not  take  time  steps  more  than  about  three 
times  larger  than  the  convective  stability  limit  due  to 
U(z) ]  . 

A  run  using  32  x 32  x 33  modes  (M=N=16,P=32) 
to  represent  each  component  of  velocity  requires  about  6s  (2.5s 

for  evaluation  of  v  x  ui)  on  the  CDC  7600  computer  per  time 
step,  including  input-output  overhead.  The  computer  time 
is  nearly  proportional  to  MNP.  A  typical  run  involves 
about  1000  time  steps.  Some  runs  to  validate  our  code 
are  reported  in  Sec.  3.  We  note  here  that  it  is  our 
(unpleasant)  experience  that  typical  transition  calculations 
require  about  10  times  as  many  time  steps  as  typical 
turbulence  decay  calculations  (see  Orszag  &  Patterson 
1972) .  The  reason  seems  to  be  the  need  to  maintain 
accurate  phase  relations  over  many  linear  oscillation  times 
in  transition  calculations  whereas  turbulent  flows  are 
nearly  critically  damped  and  evolve  quite  rapidly  in 
interesting  ways.  We  also  note  that  less  than  15  h 
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•  of  CDC  7600  time  was  spent  on  the  present  series  of 
computations  (including  a  large  number  of  runs  not 
reported  here) .  On  the  newly  introduced  CRAY-1  computer 

*  the  same  series  of  runs  requires  less  than  1.4h,  showing 
the  great  strides  that  can  now  be  made  on  significant 

fluids  problems  with  modest  computer  resources. 

t 


♦ 


* 
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I 
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3.  Two-Dimensional  Finite -Amplitude  Disturbances 

In  this  Section,  we  present  results  obtained  using 
the  computer  code  described  in  Sec.  2  for  the  evolution  of 
small  and  finite-amplitude  two-dimensional  disturbances  of  plane 
Poiseuille  and  plane  Couette  flow.  The  first  series  of 
runs  were  made  primarily  to  verify  the  accuracy  of  our  computer 
program  in  comparison  with  known  results.  Two  types  of 
tests  have  been  made:  comparison  of  the  evolution  of  small- 
amplitude  disturbances  with  their  predicted  behavior  according 
to  the  Orr-Sommerfeld  equation;  and  comparison  of  some  special 
finite-amplitude  two-dimensional  solutions  with  the  numerical 
results  of  George,  Heliums  &  Martin  (1974) . 

In  Table  1,  we  present  results  of  small -amplitude  tests 
of  the  computer  code.  The  values  of  oj  are  obtained  by 
Chebyshev-spectral  solution  of  the  Orr-Sommerfeld  equation 
and  represent  the  least  stable  eigenmode  of  the  flow  for  the 
given  a  and  8  .  The  predicted  amplitude  change  in  the 

time  interval  0<t<T  is  exp[(Im  w)T]  ,  while  the  predicted 
phase  change  in  the  same  time  interval  is  (in  radians) 

-(Re  <jj) T  .  It  seems  that  most  of  the  small  error  between 
the  predicted  and  computed  amplitude  and  phase  changes  of 
both  the  two- and  three-dimensional  modes  of  plane  Pouiseuille 
flow  and  plane  Couette  flow  is  due  to  time  differencing 
error.  Similar  tests  of  the  computer  code  on  small-amplitude 
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disturbances  have  been  made  for  Reynolds  numbers  up  to 
50,000  with  similarly  good  results. 


Our  next  verification  run  used  a  large  two-dimensional 
disturbance  studied  previously  by  George  et  al  (1974).  George 
et  al  solved  the  two-dimensional  Navier-Stokes  equations  by 
expanding  the  velocity  field  into  a  Fourier  series  in  x 
and  applying  finite-difference  methods  in  z  .  George  et  al 
did  not  use  initial  conditions  corresponding  to  eigenmodes 
of  the  Orr-Sommerfeld  equation,  but  rather  chose  the  initial 
velocity  field  to  be  that  generated  by  the  streamfunction 


<Mx,z) 


.  ,cosh  az 
'cosh  a 


cos  az  . 
cos  a 


cos  ax. 


(3.1) 


Here  a  *  2.365  in  order  to  satisfy  the  boundary  conditions, 
and  we  integrated  the  neutrally  stable  case  found  by  George 
et  al  in  which  R  =  4000,  a  =  1.05,  k  =  0.0599.  The  initial 
maximum  amplitude  of  the  perturbation  of  x-velocity  is  0.1465. 
Using  8  Fourier  inodes  in  x  and  33  Chebyshev  polynomials  in 
z  ,  the  results  of  our  computer  code  agreed  well  through  late 
times  with  those  of  George  et  al .  It  seems  that  the  reason  that 
George  et  al  were  unable  to  achieve  two-dimensional  instability 
at  Reynolds  numbers  as  low  as  those  predicted  by  the  Stuart- 
Watson  technique  is  due  to  two  causes:  (a)  the  initial 
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condition  (3.1)  is  not  a  pure  mode  of  the  Orr-Sommerfeld 
equation;  and  (b)  more  importantly,  the  choice  a  =  1.05, 
while  close  to  the  most  unstable  a  for  small-amplitude 
disturbances  is  far  from  the  most  unstable  a  for  finite- 
amplitude  disturbances  (see  below) . 

Some  characteristics  of  the  two-dimensional  finite-amplitude 
runs  discussed  below  are  given  in  Table  2.  Additional  details 
are  given  by  Kells  (1978) .  The  results  of  Herbert  (1976)  suggest 
that  the  critical  Reynolds  number  for  finite-amplitude  two-dimen¬ 
sional  instability  of  plane  Poiseuille  flow  is  roughly  2935  with 
the  unstable  mode  having  a  wavenumber  a  =  1.3231.  Runs  2-4 

present  results  of  our  numerical  simulations  of  such  a  flow. 

In  Fig.  1  we  plot  the  profile  of  the  x-velocity  component 
of  the  two  dimensional  primary  disturbance  mode  [which  depends  on 

x  like  exp(iax)]  imposed  at  t  =  0  in  Run  2.  The  initial 
disturbance  imposed  in  Runs  3  and  4  has  the  same  shape  but 
is  reduced  in  amplitude.  In  Fig.  2  we  plot  the  time  evolution 
in  Run  2  of  the  maximum  amplitude  A  of  the  x-velocity  of 
the  primary  wave  and  its  harmonic  [that  depending  on  x  like 
exp(2iax)].  After  an  initial  transient  period  the  primary  wave 
settles  down  into  a  period  of  slow  growth  suggesting  that  the 
initial  finite-amplitude  disturbance  does  not  die  out  as 
t-*00  .  However,  despite  the  growth  of  the  finite-amplitude 
disturbance,  there  is  no  sign  of  ’turbulence'  in  the  sense  that 
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the  flow  remains  ordered  and  well  behaved.  In  Fig.  3,  we  plot 
the  profile  of  the  primary  wave  after  evolution  to  t  =  120 
in  Run  2;  in  this  Figure,  the  x-velocity  component  u(z)  is 
plotted  as  a  function  of  z  at  a  point  x  at  which  its  phase 
relative  to  the  initial  perturbation  is  0  .  Comparision  of 
Figs.  1  and  3  shows  that  nonlinear  effects  tend  to  move  the 
maximum  perturbation  velocity  away  from  the  walls  at  z  =  ±1 

In  Fig.  4,  we  plot  the  unperturbed  velocity  profile  1  -  z 
and  the  mean  velocity  profile  obtained  from  Run  2  at  t  =  120. 

In  Fig.  5,  we  plot  the  curvature  of  the  unperturbed  flow 
and  the  mean  flow  at  t  =  120.  Observe  that  although  the  mean 
flow  is  changed  only  slightly  in  evolution  from  t  =  0  to 
t  =  120,  there  are  large  changes  in  the  curvature  in  the  wall 
regions . 

Run  2  has'  32  Fourier  modes  in  x  and  33  Chebyshev  modes 
in  z.  Its  numerical  accuracy  was  tested  by  making  other  runs 
with  different  values  of  M  and  P  .  A  run  with  16  Fourier 
modes  in  x  and  33  polynomials  in  z  gave  nearly  identical 
results,  while  a  run  with  8  Fourier  modes  indicated  faster  in¬ 
stability.  Increasing  the  number  of  Chebyshev  polynomials  in 
z  to  65  gave  negligible  change.  In  Run  2,  higher  harmonics 

have  Very  small  amplitudes;  the  harmonic  exp(3iax)  has  maximum 
amplitude  0.01  in  z  at  t  =  120  while  the  harmonic 
exp(7iax)  has  amplitude  0.001. 

The  results  of  Runs  3  and  4  indicate  the  effect  of  chang¬ 
ing  the  initial  amplitude  of  the  finite-amplitude  disturbance. 

In  Run  3  a  nearly  stationary  disturbance  is  achieved  in 


evolution  to  t  =  100  (see  Fig.  10) .  On  the  other  hand, 

Run  4  shows  that  there  is  a  critical  amplitude  of  about 
10%  below  which  finite-amplitude  disturbances  decay  at 
R  =  2935.  In  Fig.  6,  we  plot  the  maximum  amplitude  of  the 
x-velocity  of  the  primary  disturbance  and  its  harmonic  as 

a  function  of  time  in  Run  4.  After  an  initial  transient  period, 
both  the  primary  and  its  harmonic  settle  down  to  a  state  of 
steady  decay . 

In  Fig.  7  we  plot  the  primary  wave  and  harmonic  wave 
amplitude  vs.  time  for  Run  1  at  R  =  2500.  The  result  that 
the  disturbance  decays  as  t  increases  for  an  initial  disturb¬ 
ance  which  is  designed  to  be  close  to  the  most  unstable  finite- 
amplitude  disturbance  at  this  Reynolds  number  suggests  that  all 
two-dimensional  disturbances  decay  at  R  =  2500.  We  conclude 
that  the  threshold  for  neutral  stability  is  near  R  =  2800, 
in  rough  agreement  with  previous  theoretical  results. 

Run  5  shows  that  decreasing  a  to  a  =  1  gives  stable 
finite-amplitude  results  even  at  R  =  3500.  This  result  explains 
why  George  et  al  (1974)  were  unable  to  find  two-dimensional 
instabilities  at  Reynolds  numbers  below  3500.  George  et  al 
used  a  disturbance  wave  number  a  =  1.05  which  is  close  to 
the  value  of  a  that  gives  the  linearly  least  stable  mode 
of  plane  Poiseuille  flow  at  R  =  3500  among  all  real  a  , 
but  a  =  1.05  is  far  from  those  modes  that  give  unstable 
finite-amplitude  disturbances.  While  the  amplitudes  of  the 
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harmonics  of  the  primary  wave  disturbances  in  Runs  2  and  3  are 

quite  small,  nonlinear  effects  do  cause  considerable  distortion 

of  the  primary  wave  itself  (see  Figs.  1  and  3) ,  giving  a  large 

effect  on  the  least  stable  a  . 

In  summary,  while  there  are  finite-amplitude  two-dimensional 

instabilities  of  plane  Poiseuille  flow  at  subcritical  Reynolds 

numbers,  there  are  apparently  no  explosive  instabilities  that 

can  generate  small-scale  random  flow  structures  when  only 

two-dimensional  interactions  are  allowed. 

Run  6  illustrates  the  effect  of  finite-amplitude  two- 

dimensional  disturbances  on  plane  Couette  flow.  In  Fig.  8 

we  plot  the  profile  of  the  x-velocity  component  of  the  initial 

disturbance  in  Run  6.  This  initial  disturbance  is  constructed 

from  a  solution  of  the  Orr-Sommerfeld  equation  as  follows. 

For  the  given  a  =  2,  8=0,  and  R  =  5000,  the  least  stable 

A 

eigenmode  u  (z)  of  plane  Couette  flow  is  asymmetric  about  z=0. 

It  is  easy  to  verify  from  the  Orr-Sommerfeld  equation  (1.6) 

A 

that  if  u(z)  corresponds  to  the  wave  vector- frequency 
a, 8, m,  then  u  ( —  z )  is  the  eigenfunction  associated  with 
a, 8,  -Re (w) +ilm(w) .  Note  that  the  phase  velocity  of  this 
complex- congugate  reflected  mode  is  the  negative  of  the  phase 
velocity  of  the  original  mode.  The  initial  conditions  for 
Run  6  are  chosen  to  be  the  (symmetrized)  conditions  in  which 
the  x-velocity  is  given  by 

u(x,z,0)  =  ^Re  [u(  z) +u*  (-z)  ] eiax.  (3.2) 

In  Fig.  9  we  plot  the  time  evolution  of  the  maximum  amplitude 

of  the  primary  disturbance  and  its  harmonic  as  functions  of  time 

for  Run  6.  The  decay  is  quite  rapid,  in  contrast  with  the 
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predictions  of  Ellingsen,  Gjevik  &  Palm  (1970)  and  Coffee 
(1977)  whose  theoretical  calculations  seem  to  predict  that 
Run  6  should  lead  to  a  finite-amplitude  stationary  disturbance. 

Other  runs  for  plane  Couette  flow  give  '.o  indication  as  yet 
of  any  region  of  two-dimensional  instability  [although  our 
computer  code  is  not  capable  of  resolving  the  case  suggested 
by  Kuwabara  (1967) ] .  It  is  not  likely  that  our  decay  results 
for  plane  Couette  flow  are  due  to  the  low  8  *33  resolution 
used  in  Run  6.  The  33  Chebyshev  polynomials  used  in  z  are 
more  than  adequate  to  resolve  the  vertical  structure.  The 
8  Fourier  modes  used  in  x  also  seem  to  be  adequate;  at 
t  =  90,  the  primary  wave  has  maximum  amplitude  0.0053,  its 
harmonic  has  maximum  amplitude  0.00034,  while  its  second 
harmonic  has  amplitude  0.000060.  In  any  case,  our  experience 
with  increasing  the  x-resolution  is  that  higher  resolution 
gives  more  stable  results  than  low  resolution  (see  Sec.  4) . 

The  high  frequency  oscillation  exhibited  by  the  primary 
wave  in  Fig.  9  is  due  to  the  fact  that  the  maximum  amplitude 
of  u(z)  appears  alternately  near  z  =  ±1;  this  effect  is 
a  consequence  of  the  symmetrized  mode  (3.2)  in  which  the  two 
components  have  equal  and  opposite  phase  velocity.  Unsymmetrized 
initial  conditions  do  not  lead  to  these  rapid  oscillations (see  Fig.  37) . 

Our  results  suggests  a  serious  discrepancy  between  the 
Stuart-Watson-Eckhaus  theory  predictions  of  two-dimensional 
finite-amplitude  instability  of  plane  Couette  flow  and  direct 
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numerical  calculations  that  indicate  decay.  One  possible  explanation 

is  that  we  have  not  sought  solutions  with  large  enough  a. 

According  to  Davey  &  Nguyen  (1971)  ,  the  finite-amplitude  disturbances 

1/2 

requiring  the  least  energy  to  excite  have  a  -  0.1 3R  '  or 
a~9  for  R  =  5000.  However,  the  energy  required  to  excite  a 
finite-amplitude  mode  at  this  large  Reynolds  number  is  a  fairly 
flat  function  of  a  and  our  inability  to  find  instability  with 
a=2  is,  in  our  view,  a  serious  criticism  of  the  theory.  We 
are  now  studying  solutions  with  larger  a  using  a  higher  resolution 
computer  code  in  order  to  resolve  the  discrepancy  between  non¬ 
linear  stability  theory  and  numerical  experiments.  However,  the 
main  point  of  the  present  paper  is  the  strength  of  three- 
dimensional  effects;  the  possibility  that  two-dimensional  inter¬ 
actions  may  be  slightly  stronger  for  much  larger  a  than 
considered  here  is,  we  believe,  a  secondary  issue. 

In  summary,  our  direct  calculations  of  the  Navier-Stokes 
equations  are  in  reasonably  good  agreement  with  nonlinear  stab¬ 
ility  theory  calculations  of  plane  Poiseuille  flow.  While  no 
transition  to  turbulence  is  either  observed  or  predicted  due  to 
two-dimensional  disturbances,  there  are  finite-amplitude  motions 
that  do  not  decay  for  Reynolds  numbers  larger  than  about  2800. 

Serious  discrepancies  do  now  exist  between  nonlinear  stability 
theory  and  direct  calculations  of  two-dimensional  finite-amplitude 
effects  in  plane  Couette  flow. 
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Three-Dimensional  Finite-Amplitude  Disturbances 


In  this  Section,  we  present  results  obtained  using  the 

computer  code  described  in  Sec.  2  for  the  evolution  of  finite- 

amplitude  three-dimensional  disturbances  of  plane  Poiseuille 

and  plane  Couette  flow.  Some  characteristics  of  the  runs 

discussed  in  this  Section  are  given  in  Table  3  for  the  plane 

Poiseuille  runs  and  in  Table  4  for  the  plane  Couette  flow  runs. 

Run  3A  is  identical  to  Run  3  discussed  in  Sec.  3  except 

for  one  interesting  three-dimensional  feature.  The  two- 

dimensional  Run  3  was  made  with  the  three-dimensional  computer 

code  using  8  x  4  x  33  resolution  and  setting  all  flow  components 

to  be  independent  of  y  at  t  =  0.  Since  we  did  not  reset  the 

flow  to  be  independent  of  y  at  later  times,  round-off  error 

-14 

induced  a  small  (order  10  )  three-dimensionality  on  the 

flow.  This  three-dimensional  component  was  rapidly  filtered 

into  an  unstable  three-dimensional  mode  that  grew1  rapidly  on 

the  finite-amplitude  two-dimensional  background.  A  plot  of  the 

maximum  amplitude  of  the  primary  two-dimensional  disturbance, 

its  two-dimensional  harmonic,  and  the  fundamental  oblique  wave 

(with  wave  numbers  a  =  1.3231,  0=1)  is  given  in  Fig.  10. 

Note  that  the  amplitude  of  the  oblique  component  is  multiplied 
0 

by  ID  .  The  results  of  Run  3A  show  that  the  evolved  two- 
dimensional  f inite-amj litude  state  is  quite  susceptible  to  three- 
dimensional  instabilities. 

Run  3A  suggests  a  simple  mechanism  for  transition.  Whereas 
two-  and  three-dimensional  linear  instabilities  of  plane 


Poiseuille  flow  are  quite  mild  and  occur  only  for  large 
Reynolds  numbers  (R^_577 2),  finite-amplitude  two-dimensional 
neutrally  stable  states  of  plane  Poiseuille  flow  seem  to 
be  explosively  unstable  to  three  dimensional  perturbations 
at  subcritical  Reynolds  numbers. 

In  Runs  7-16  we  investigate  the  finite-amplitude  aspects 
of  three-dimensional  instability  of  plane  Pouiseuille  flow. 

All  of  these  runs  (except  Run  16)  are  made  by  choosing  the 
initial  condition  to  consist  of  a  superposition  of  four 
components : 


V  =  (U(z)  ,0,0)+Retv2D(z)eiaX]+Re[v+^3D(z)eiaX+lBY) 


_  .  ,  .  iax-igy. 

+Re[v_  3D(z)e  1  ] 


(4.1) 


Here  U(z)  =  1-z2  is  the  unperturbed  plane  Pouiseuille  flow, 

v  (z)  is  a  two-dimensional  eigenf unction  of  the  Orr-Sommerfeld 
-2D' 

equation,  and  v+  3D(z)  is  a  three-dimensional  eigenfunction 

of  the  Orr-Sommerfeld  equation  with  spanwise  wave  number 

±6  .  The  form  of  this  initial  flow  field  is  chosen  to  corresp 

ond  to  that  suggested  by  the  experiments  of  Klebanoii,  Tidstrom  S. 

Sargent  (1962)  and  the  theory  of  Benney  &  Lin  (1960) .  The  oblique 

three-dimensional  .disturbance  is  always  formed  with  symmetric 

±(J  components,  so  the  initial  conditions  give  a  standing 

wave  disturbance  with  the  x-axis  an  axis  of  symmetry.  For 

all  runs  (except  Runs  12  and  13)  the  Orr-Sommerfeld  eigenfunctions 
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~2D  anc^  ^3D  are  c^osen  to  be  the  least  stable  eigenfunctions 
for  the  given  values  of  a  ,3  .  Because  of  the  computational 
cost  of  the  three-dimensional  runs,  we  have  not  been  able  to 
make  a  complete  survey  of  the  effects  of  disturbance  wave 
length  and  amplitude  on  the  strength  of  the  transition  process. 
However,  there  are  indications  on  the  basis  of  the  results 
reported  below  that  the  chosen  values  of  a  and  0  give 
nearly  the  most  strongly  unstable  results. 

In  Figs.  11  and  12  we  plot  the  profiles  of  the  two  dimensional 
disturbance  u2D(z)  an<^  three-dimensional  disturbance  u+3D(z), 

respectively,  applied  initially  in  Runs  7  and  8.  In  Figs. 

13  and  14,  we  plot  the  maximum  amplitude  of  the  two-dimensional 
primary  disturbance,  its  two-dimensional  harmonic  [depending 
on  x  like  exp(2iax)]  ,  and  the  primary  oblique  wave  [depending 
on  x  and  y  like  exp(iax+if3y)  ]  for  Runs  7  and  8,  respectively. 

The  only  difference  between  Runs  7  and  8  is  that  Run  7  has  twice 
the  x-y  resolution.  It  is  apparent  from  Figs.  13  and  14  that 

the  results  are  nearly  identical  through  t  =  35,  but  the  results 
later  diverge.  In  Run  7,  the  flow  seems  to  'break  down'  at 
t  ~  40  ,  while  in  Run  8  the  'breakdown'  seems  to  occur  for 
t  =  37.  (Part  of  this  difference  may  be  due  to  time  truncation 
error  since  Run  7  uses  a  time  step  of  0.05  while  Run  8 
uses  a  time  step  of  0.1.)  Another  run  made  with  the  initial 
conditions  of  Runs  7  and  8  but  increasing  the  z-resolution 
to  65  Chebyshev  polynomials  showed  no  change  from  the  results 
with  33  Chebyshev  polynomials  until  well  into  the  breakdown 


region. 


In  Tables  5  and  6,  we  give  values  of  the  maximum  amplitudes 
1  n  z  t^ie  x-velocity  Fourier  components  u(m,n,z,t),  given 
in  terms  of  the  spectral  components  u(m,n,p,t)  in  (2.3)  by 

P 

u(m,n,z,t)  =  l  u (m, n, p, t)  T  (z) 

p=0  p 

While  the  convergence  of  the  expansion  (2.3)  evidently  does 
deteriorate  somwehat  with  time,  especially  near  t  ~ 40 
when  the  flow  breaks  down,  the  decrease  of  the  maximum  amplitude 
with  increasing  m  and  n  is  sufficiently  rapid  for  us  to 
confidently  assert  that  the  results  of  Run  7  are  real  and  not 
due  to  numerical  instability.  Also,  while  the  effects  of  truncation 
of  the  Fourier  expansion  (2.3)  in  Run  8  are  even  more  pronounced 
than  for  Run  7,  the  lower  resolution  results  of  Run  7,  are  a  faithful 
predictor  of  flow  breakdown  and  accurately  represent  the  flow 
until  just  before  breakdown  occurs. 

We  emphasize  the  qualitative  change  in  the  behavior  of 
the  flow  brought  about  by  three  dimensionality.  While  the 
finite  amplitude  two-dimensional  results  presented  in  Sec.  3 
show  the  resulting  flow  to  be  smooth,  periodic  and  with  no 
appreciable  small-scale  excitation,  quite  the  opposite  is  true 
in  the  three-dimensional  flows  considered  here.  In  three 
dimensions,  the  flow  develops  into  a  state  of  apparent  random¬ 
ness  with  appreciable  small-scale  excitation.  It  is  this  latter 
state  that  we  claim  is  representative  of  the  transition  to 
turbulence . 

While  we  cannot  claim  that  our  results  for  Run  7  are  in  detail 
accurate  for  t  >50  when  the  flow  is  seemingly  random  and 
"turbulent",  the  relative  insensitivity  to  changes  in  resolution 
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suggest  strongly  that  our  simulations  accurately  portray  the 
breakdown  of  the  laminar  flow  and  its  transition  to  turbulence. 
[While  the  flow  after  breakdown  is  almost  certainly  not  being 
treated  in  detail  accurately,  there  is  some  basis  for  the 


assertion  that  the  statistical  properties  of  the  resulting  flow 
are  in  reasonable  agreement  with  real  flows.  Some  further 
elaboration  of  this  point  is  given  below.  However,  we  do  not 
stress  the  turbulence  aspects  of  the  developed  flow  in  this 
paper  as  numerical  simulations  of  turbulence  are  best  done 
slightly  differently.] 

The  apparent  discontinuities  in  the  time-history  curves 
plotted  in  Figs.  13,  14,  and  later  arise  because  the  maximum 
amplitude  of  each  of  the  flow  components  may  come  from  different 
local  maxima  and,  hence,  quite  distinct  z  locations  at  different 
times . 


In  Fig.  15,  we  plot  x-y  averages  of  the  x-component 
of  the  velocity  field  in  Run  8  at  t  =  30-75(15): 
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The  results  for  Run  7  at  t  =  30  and  t  =  45  (the  only  times 
at  which  corresponding  results  are  available)  are  indistinguishable 
from  the  plotted  curves.  The  mean  velocity  profiles  plotted  in 
Fig.  15  correspond  to  profiles  in  the  nonlinear  laminar. 
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early  transition,  late  transition,  and  turbulent  flow  regimes. 

The  mean  velocity  profile  at  t  =  75  is  strikingly  similar  to 
mean  velocity  profiles  observed  experimentally  in  turbulent 
channel  flow  (Laufer  1951,  Comte-Bellot  1965) .  In  a  continuation 
of  Run  8  to  t  >75,  no  change  in  the  mean  velocity  profile  is 
found. 

While  there  are  no  strong  inflections  appearing  in  the 
mean  velocity  profiles  plotted  in  Fig.  15,  instantaneous 
velocity  profiles  plotted  in  Figs.  16  and  17  for  Run  7  at 
t  =  30  and  t  =  45,  respectively,  show  strong  inflections. 

These  local  velocity  profiles  are  strongly  unstable  and, 
presumably,  their  instability  leads  to  the  random  behavior 
observed  at  later  times. 

Contour  plots  of  the  x-velocity  component  for  Run  7 
in  the  y-z  plane  at  x  =  0  are  given  in  Fig.  18  for 
t  =  0-45(15)  .  At  the  spanwise-y  location  of  the  center  of 
the  'cat's-eyes'  observed  in  Fig.  lEf,  the  velocity 
profile  is  much  fuller  than  the  undisturbed  plane  Poiseuille 
profile  and  the  corresponding  wall  shear  is  much  greater  (see 
Figs.  16  and  17  )  .  In  the  region  between  the  'cat's-eyes', 
the  local  velocity  profile  is  highly  inflectional  (as  shown 
in  Figs.  16  and  17)  . 

Contour  plots  of  the  longitudinal  vorticity  component 
for  Run  7  in  the  y-z  plane  at  x  =  0  are  given  in  Fig.  19  for 
t  =  0-45(15).  Several  features  of  these  contour  plots  should 
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be  observed.  First,  there  is  considerable  enhancement  of 
longitudinal  vorticity  during  time  evolution  of  this  flow. 

In  Table  5,  we  present  data  on  the  maximum  amplitude  of  the 
transverse  waves  [m=0  in  (2.3)]  as  a  function  of  time  in  Run  7. 

Through  t  =  35,  the  data  in  Table  5  agrees  to  within  a  relative 
error  of  less  than  2%  with  the  corresponding  data  for  Run  8. 

The  primary  transverse  mode  seems  to  achieve  rapidly  a  nearly 
steady  value,  while  the  flow  breakdown  is  signaled  by  rapid 
growth  of  the  n=2  harmonic  component.  Second,  the  small- 
scale  structure  evident  in  Fig.  19  at  t  =  45  indicates  the 
breakdown  of  the  flow. 

In  Fig.  20,  contour  plots  are  given  of  the  y-  and  z- 
components  of  the  velocity  and  vorticity  for  Run  7  in  the 
y-z  plane  at  x  =  0  at  t  =  30.  Examination  of  these  contour 
plots  (and  others  made  at  different  locations)  show  several 
interesting  features.  First,  there  is  evidently  a  correlation 
between  the  x-  and  z-velocity  components  representing  enhanced 
transfer  of  momentum  towards  the  walls.  Also,  the  contours  of 
u  in  the  x-y  plane  given  in  Fig.  21  sho w  gusts  of  high- 
velocity  fluid  moving  through  the  center  of  an  otherwise 
slow-moving  stream;  the  leading  edges  of  a  gust  seem  to  be 
sharper  than  the  trailing  edges,  implying  negative  skewness 
in  velocity  derivatives. 

In  Fig.  22,  we  plot  the  profile  of  the  two- 
dimensional  primary  wave  at  t  =  30  and  t  =  45  ,  respectively, 
for  Run  7.  In  Fig.  23  ,  we  present  a  similar  plot  of  the 
two-dimensional  harmonic  components  at  t  =  45  in  Run  7, 
while  in  Fig.  24  we  plot  the  profile  of  the  primary  oblique 
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wave  (a  =  l,g  =  1)  at  t  -  45  in  Run  7.  Figs.  22-24  show 
that,  through  the  period  of  initial  breakdown,  the  largest 

fluctuations  occur  in  the  neighborhood  of  the  walls.  At  later 
times,  the  harmonic  components  develop  slightly  stronger  activity 
near  the  center  of  the  channel.  The  spatial  structure  and  amplitude 
of  these  fluctuations  is  not  inconsistent  with  experimental 
observations  (Laufer  1951  ,  Comte-Bellot  1965)  . 

Run  9  tests  the  effect  of  decreasing  the  initial  amplitude 
of  the  two-dimensional  and  three-dimensional  initial  disturbances 
in  Runs  7  and  8  by  25%.  In  Fig.  25,  the  maximum  amplitude  of  the 
primary  two-dimensional  wave,  its  harmonic,  and  the  primary  three- 
dimensional  wave  for  Run  9  are  plotted  as  functions  of  t  .  The 
instability  in  Run  9  is  considerably  weaker  than  in  Run  8.  In 
other  runs,  we  have  found  that  keeping  the  amplitude  of  the 
two-dimensional  disturbance  in  Run  8  at  0.11,  but  decreasing 
the  amplitude  of  the  three-dimensional  disturbance  to  0.025 
stabilizes  the  flow.  Similarly,  another  run  decreasing  the 
amplitude  of  the  two-dimensional  wave  to  0.05,  but  keeping  the 
amplitude  of  the  three-dimensonal  wave  disturbance  at  0.05  as 
in  Run  8,  stabilized  the  flow.  At  R  =  1250,  there  seems  to 
be  a  critical  amplitude  of  about  0.08  for  the  two-dimensional 
component  and  about  0.04  for  the  three-dimensional  component. 

Runs  10-13  were  designed  to  test  the  sensitivity  of  the 
flow  to  various  spanwise  wave  numbers.  The  amplitude  of  the 
initial  disturbances  in  these  runs  is  chosen  to  be  close  to  that 
of  Runs  7  and  8.  Run  10  uses  a  primary  spanwise  wave  number 
8  =  0.25.  In  this  case,  the  three-dimensional  disturbance  is 
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nearly  two-dimensional  and  is  barely  distinguishable  from  the 

primary  two-dimensional  disturbance.  The  resulting  flow  is 
so  nearly  two-dimensional  that  the  disturbances  die  away  quickly. 

Run  11  uses  a  primary  spanwise  wave  number  3  =  0.5  ;  in  this 
case  transition  was  observed  but  only  at  the  late  time  t  =  120. 
(Even  though  Runs  10-13  use  the  low  spatial  resolution 
8  *  8  *  33  ,  we  believe  that  the  prediction  of  transition  in 
Runs  11  and  12  is  justified  because  the  three-dimensional  distur¬ 
bance  remains  at  large  amplitude  during  the  laminar  phase.  The 
resolution  problems  discussed  below  with  regard  to  Runs  14  and  15 
seem  to  apply  only  to  marginally  unstable  flows  •) 

Other  runs  made  with  the  two-  and  three-dimensional 
disturbances  detuned  from  the  least  stable  modes  of  the  Orr- 
Sommer feld  equation [used  in  (4.1)]  gave  significantly  less 
tendency  to  undergo  transition  than  the  runs  discussed  here. 

Runs  12  and  13  introduce  a  new  feature  into  the  calculation. 
With  6=2  (Run  12)  and  6=4  (Run  13) ,  the  least  stable  three- 
dimensional  mode  of  the  Orr-Sommerf eld  equation  is  no  longer  concent¬ 
rated  near  the  wall  like  the  two-and  three-dimensional  modes  plotted 

in  Figs.  11  and  12.  The  least  stable  three-dimensional  mode  with 
6=4  at  R  =  1250  is  plotted  in  Fig.  26;  it  is  evident  that 
this  mode  is  concentrated  near  the  center  of  the  channel.  Runs 
made  using  this  three-dimensional  center  mode  together  with 
the  two-dimensional  disturbance  mode  plotted  in  Fig.  11  invariably 
decay,  evidently  because  there  is  not  enough  interaction  between 
the  center  mode  and  the  wall  mode  disturbances.  However,  transition 
to  turbulence  with  these  larger  spanwise  wave  numbers  6  can  be 
achieved  by  using  a  three-dimensional  wall  mode  disturbance. 
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The  least  stable  three-dimensional  wall  mode  disturbance  with 
8=4  at  R  =  1250  is  plotted  in  Fig.  27.  Runs  12  and  13  use 
these  least  stable  wall  modes.  Run  12  is  strongly  unstable 
and  leads  quickly  to  transition.  On  the  other  hand,  Run  13  is 

stable  and  does  not  undergo  transition,  even  though  the 

wall  mode  disturbance  is  used.  Evidently  there  is  a  spanwise 

wave  number  selection  mechanism  involved  in  the  transition 

process.  It  seems  that  the  most  dangerous  three-dimensional  disturbances 

have  wavefronts  at  an  angle  of  about  45°-60°  to  the  mean  flow. 

Runs  14  and  15  illustrate  a  possible  pitfall  of  numerical 
analysis  of  these  transition  problems.  The  only  difference  between 
these  two  runs  is  the  horizontal  resolution:  in  Run  14,  only 
8x8  Fourier  modes  are  used  to  resolve  the  x  and  y  directions; 
in  Run  15,  16  x  16  Fourier  modes  are  used.  Both  runs  are 

made  at  R  =  750.  In  Fig.  28,  we  plot  the  maximum  amplitudes 
of  the  primary  two-dimensional  wave,  its  harmonic,  and  the 
primary  three-dimensional  wave  for  Run  14.  The  disturbances 
decay  until  about  t  =  100  and  then  abruptly  erupt  into 
turbulence.  A  similar  plot  for  Run  15  is  given  in  Fig.  29. 

In  this  case,  the  disturbances  undergo  only  a  smooth  continuous 
decay  and  no  'transition'  is  observed.  We  have  also  made  a  run 
in  which  the  8  x  8  x  33  resolution  of  Run  14  is  increased  to 
16x16x33  at  t  =  90  in  Run  14;  in  this  case,  as  in  Run  15, 
no  transition  is  observed.  Evidently  low-resolution  can  give 
premature  predictions  of  transition.  Similar  effects  of  low 
resolution  have  been  observed  by  McLaughlin  &  Orszag  (1979)  in 
calculations  of  transition  in  three-dimensional  Benard 
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convection.  It  is  essential  that  transition  calculations  be 


done  with  utmost  care,  in  order  to  avoid  these  spurious  predictions 
of  instability  and  breakdown  in  three-dimensional  calculations. 

Our  spectral  calculations  do  provide  a  'bootstrap'  procedure  to 
test  internal  accuracy  of  the  simulations;  if  the  spectrum 
obtained  by  Run  14  is  tested  for  convergence  in  the  same  way 
as  the  spectrum  of  Run  7  is  tested  by  the  results  given  in 
Tables  5  and  6,  it  is  found  that  the  8><8  truncation  in 
x  and  y  in  Run  14  has  a  large  effect. 

We  are  confident  that  our  transition  predictions  cited  earlier 
in  this  Section, especially  for  Run  7,  are  not  affected  by  these 
resolution  problems  of  Run  14. 

Our  final  plane  Poiseuille  flow  run  is  Run  16.  Here  the 
initial  conditions  consist  of  the  output  from  Run  8.  These 
"turbulent"  initial  conditions  are  then  run  at  the  lower 
Reynolds  number  R  =  500,  in  order  to  test  whether  an  initial 
field  of  turbulence  can  persist  at  this  lower  Reynolds  number. 

In  Fig.  30,  we  plot  the  evolution  of  the  maximum  amplitude 
of  the  two-dimensional  primary  wave,  its  harmonic,  and  the  primary 
three-dimensional  wave.  While  the  results  are  not  conclusive, 
it  seems  that  the  field  of  turbulence  is  slowly  decaying  and 
that  R  =  500  is  not  sufficiently  high  to  sustain  turbulence. 

Runs  17-20  investigate  the  effect  of  three-dimensional  finite- 
amplitude  disturbances  on  plane  Couette  flow.  While  plane  Couette 
flow  is  much  more  stable  than  plane  Poiseuille  flow  for  small- 
amplitude  disturbances  and  two-dimensional  finite- 
amplitude  disturbances  (see  Sec.  3),  the  effects 
of  finite-amplitude  three-dimensional  disturbances  at  modest 


Reynolds  numbers  are  just  as  dramatic  in  plane  Couette  flow  as 
in  plane  Poiseuille  flow. 


In  Figs.  31  and  32,  we  plot  the  profiles  of  the  x-velocity 
of  the  two-  and  three-dimensional  primary  disturbances,  respectively, 
applied  initially  in  Runs  17  and  18  at  R  =  1250.  These  profiles 
are  obtained  from  the  corresponding  least  stable  eigenfunctions 
of  the  Orr-Sommerfeld  equation  for  plane  Couette  flow  by  forming 
the  symmetric  combination  (3.2)  [with  an  extra  factor  exp(igy) 
for  the  three-dimensional  mode] .  (The  asymmetry  observed  in 
Fig.  31  is  due  to  the  complex  phase  of  the  mode  at  the  particular 
value  of  x  at  which  the  plot  was  made.)  In  Fig.  33,  we  plot 
the  maximum  amplitude  of  the  two-dimensional  primary  disturbance, 
its  harmonic,  and  the  primary  three-dimensional  disturbance  vs 

t  for  Run  17.  The  flow  breaks  down  to  turbulence  near  t  =  45  . 

A  similar  plot  for  the  lower  resolution  Run  18  is  given  in  Fig. 

34.  There  is  good  agreement  between  the  results  plotted  in  Figs. 

33  and  34  until  beyond  the  breakdown  of  the  laminar  flow.  We 
conclude  that  this  flow  does  undergo  transition  to  turbulence. 

In  Fig.  35,  we  plot  the  mean-velocity  profile  u(z) 
for  Run  17  at  t  =  60  .  It  is  apparent  that  the  mean- 
velocity  profile  is  tending  toward  the  characteristic 
S-shape  expected  in  turbulent  Couette  flow. 

The  effect  of  the  symmetrized  initial  condition  (3.2)  on  the 
evolution  of  the  flow  is  illustrated  by  Runs  19  and  20  made  at 
R  =  1000.  The  maximum  amplitude  plot  for  Run  19  plotted  in  Fig. 

36  suggests  possible  transition  near  t  =  75  ;  we  hesitate  to 
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claim  that  this  transition  is  real  because  of  possible  resolution 
limitations.  The  maximum  mode  amplitudes  for  Run  19  listed  in 
Table  7  show  that  the  transverse  (m  =  0)  modes  are  likely  to 
be  inadequately  resolved  for  t  >  60.  An  interesting  feature 
of  Figs.  33,  34,  and  36  for  the  symmetrized  initial  disturbances 
is  the  high  frequency  oscillation  of  the  maximum  amplitude  of 
the  two-dimensional  primary  wave  and  its  harnomic.  Run  20  is 
made  using  unsymmetrized  Orr-Sommerfeld  eigenfunctions  as 
initial  conditions.  The  maximum  amplitude  plot  for  this  run,  made 
at  R  =  1000,  is  given  in  Fig.  37.  Two  features  are  notev/orthy. 
First,  in  contrast  with  Run  19,  the  transition  in  this  case  can 
be  much  more  confidently  asserted,  because  the  resolution  limit¬ 
ations  do  not  appear  to  be  severe  until  after  breakdown  occurs. 
Second,  the  high  frequency  oscillations  in  the  two-dimensional 
disturbances  and  its  harmonic  have  disappeared. 

As  for  plane  Poiseuille  flows,  the  plane  Couette  flow 
runs  are  characterized  by  very  rapid  generation  of  the  transverse 
Fourier  components  m  =  0  in  (2.3)  [see  Table  7  for  the  results 
of  Run  19].  In  these  flows,  it  seems  that  small-scale  structures 
are  generated  by  the  strong  instability  of  the  flows  resulting 
from  superposition  of  the  longitudinal  vortices,  represented 
by  the  transverse  (  m  =  0)  Fourier  components,  on  the  basic 
laminar  flows  (1.3-4). 

We  have  encountered  significant  difficulty  in  extending 
our  plane  Couette  flow  runs  to  later  times  than  initial  breakdown. 
Evidently,  the  turbulence  that  develops  is  of  a  particularly 
severe  kind  that  is  inadequately  resolved  using  the  current 
codes. 

In  conclusion,  it  seems  that  plane  Couette  flow  undergoes 
transition  at  Reynolds  numbers  at  least  as  low  as  those  for  which 
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plane  Poiseuille  flow  undergoes  transition.  Three-dimensional 
effects  are  crucial  in  establishing  breakdown  at  Reynolds  numbers 
of  order  1000. 

5.  Conclusions 

The  results  presented  in  Sects.  3  and  4  show  the  central 
role  played  by  the  interaction  of  two-and  three-dimensional 
finite-amplitdue  disturbances  in  the  breakdown  of  plane  Poiseuille 
and  plane  Couette  flows.  The  basic  character  of  this  interaction 
is  qualitatively  consistent  with  the  theory  developed  by  Bennev  & 
Lin(1960)  to  explain  the  experiments  of  Klebanoff,  Tidstrom  & 

Sargent  (1962).  However,  as  has  been  emphasized  by 
Stuart  (1961),  the  Benney-Lin  theory  can  be  at 

best  qualitatively  correct.  The  theory  assumes  that  the  phase 
velocities  of  the  interacting  two-  and  three-dimensional  waves 
are  identical,  which  is  not  correct  (see  Tables  3  and  4) . 

Our  calculations  show  that  three-dimensional  finite- 
amplitude  effects  produce  strong  inflectional  velocity  profiles 
that  eventually  break  down  to  turbulence.  In  plane  Poiseuille 
flow,  these  three-dimensional  effects  due  to  initial  disturbances 
with  amplitudes  of  5-10%  of  the  mean  flow  explain  the  experimentally 
observed  transitions  at  Reynolds  numbers  of  order  1000,  whereas 
arbitrarily  large  two-dimensional  finite-amplitude  disturbances 
seem  powerless  at  Reynolds  numbers  much  below  3000.  The  most 
dangerous  three-dimensional  interactions  seem  to  be  between 
oblique  Orr-Sommerf eld  modes  propagating  at  about  45°  to  the 
unperturbed  flow  with  a  wavelength  about  3  times  larger  than  the 
channel  depth.  In  plane  Couette  flow,  our  numerical  results 
suggest  that  three-dimensional  effects  due  to  initial  disturbances 
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with  amplitudes  of  order  5-10%  can  drive  transition  at  Reynolds 
numbers  of  order  1000,  while  two-dimensional  effects  do  not  seem 
strong  even  at  Reynolds  numbers  an  order  of  magnitude  larger. 

The  numerical  results  also  suggest  that  turbulence  can  be 
sustained  in  these  planar  shear  flows  at  somewhat  lower  Reynolds 
numbers  but  not  as  low  as  500. 

One  disturbing  feature  is  the  high  resolution  in  both 
space  and  time  that  seems  to  be  necessary  to  compute  these  trans 
ition  flows.  To  compute  transition  accurately,  it  is  necessary 
to  calculate  relatively  weak  interactions  over  many  linear 
oscillation  periods.  While  turbulent  flows  do  require  high 
spatial  resolution,  they  also  evolve  quickly  so  that  the  total 
computer  time  may  be  less  than  for  a  transition  calculation. 
Also,  in  turbulence  calculations,  only  statistical  averages 
need  be  determined  accurately  and  practice  has  shown  that 
accurate  statistical  results  can  be  obtained  with  relatively 
low  resolution. 

It  is  interesting  that  we  have  found,  in  contrast  to 
some  previous  investigators,  that  the  accuracy  requirements 
of  transition  calculations  are  more  severe  in  the  horizontal 
x-  and  y-directions  than  in  the  z-direction  normal  to  the 
walls.  The  Chebyshev  expansions  used  in  the  z-direction 
have  extraordinarily  good  resolution  near  the  walls.  Our 
result  that  low  horizontal  resolution  can  give  spurious 
predictions  of  transition  must  be  considered  carefully  in 
future  work  on  these  problems.  Low  horizontal  resolution 
prevents  the  excitation  of  small  scale  motions  that  can  act 
as  an  'eddy'  viscosity  that  damps  out  instabilities. 
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Table  1.  Behavior  of  small-amplitude  disturbances  computed 
numerically  compared  with  behavior  predicted  by 
Orr-Sommerfeld  equation. 


i 

Plane  Poiseuille  flow 

Plane  Couette  flow 

Two- 

dimensional 

disturbance 

Three- 

dimensional 

disturbance 

Two- 

dimensional 

disturbance 

Three- 

dimensional 

disturbance 

Reynolds  number  R 

1500 

1500 

5000 

5000 

x-wave number  a 

1 

1 

0.5 

0.5 

y-wavenumber  8 

0 

1 

0 

0.5 

f.e  w 

0.3262988 

0.4012928 

0.3511072 

0.3517084 

Im  w 

-0.0282057 

-0.0282305 

-0.0413797 

-0.0418436 

Initial  amplitude 
^ (x-velocity) 

2.198  xlO"5 

4 . 9 64  x  10-5 

1.209  x 10"4 

6.338  x 10"4 

Spatial  resolution 
(2M)  MP+1) 

8  x  33 

8  x  33 

8  x  65 

CO 

> 

cr> 

Time  step  At 

0.1 

0.1 

0.1 

0.1 

i’inal  time  T 

10 

10 

6 

6 

Computed  amplitude 
decay  (Kt<T 

0.753519 

0.756756 

0.781904 

0.778925 

Predicted  amplitude 
'decay  exp[(Im  u>)T] 

0.754231 

0.754044 

0.780143 

0.777975 

Computed  phase 
change  (radians) 

0<t<T 

-3.25644 

-4.00617 

-2.10548 

-2.10545 

t 

Vredicted  phase 
change  “(Re  <*>)T 

« 

1 

-3.26299 

-4.01293 

-2.10664 

-2.11025 

5G 


Table  2.  Finite-amplitude  two-dimensional  disturbance 
characteristics  of  runs  discussed  in  Sec.  3. 


Table  4.  Finite-amplitude  three-dimensional  disturbance  characteristics  of  plane 

Couette  flow  runs  discussed  in  Sec.  4. 
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Table  6. 


Maximum  Amplitude  in  z  of  the  Fourier  components  of  the  x-velocity 
u(m/n/z/t)  in  the  mixed  spectral  (x-y) -physical  (z)  space 
representation  for  Run  7. 


Time  (t) 


Component  (m,n) 

20 

30 

40 

50 

Two  Dimensional 

Harmonics 

(1,0) 

5 . 69 (-2 ) 

4.11  (-2) 

2.70  (-2) 

1 . 78  (-2) 

(2,0) 

4.50  (-3) 

4 . 45  (-3) 

3. 31 (-3) 

7.66 (-3) 

(3,0) 

1.45  (-3) 

8.31  (-4) 

3.05  (-4) 

5 . 33  (-3) 

(4,0) 

4.99  (-4) 

9.63  (-5) 

2.16  (-4 ) 

3 . 75  (-3) 

(6,0) 

6. 08  (-5) 

3.02  (-5) 

3. 50 (-5) 

2 . 26 (-3 ) 

(8,0) 

1 .10  (-5) 

1.03  (-5) 

1 . 28 (-5) 

5.92 (-4 ) 

Oblique 

Harmonics 

(1,1) 

7.61 (-2) 

8.15  (-2) 

6.14  (-2) 

2.31(-2) 

(2,1) 

6 . 75  (-3) 

2.47  (-3) 

8.71  (-4) 

9 . 4  6  (-  3) 

(4,1) 

5 . 98 (-4 ) 

2 . 80  (-4 ) 

2 . 32 (-4 ) 

1 . 89 (-3) 

(6,1) 

6.35  (-5) 

4.34  (-5) 

3 . 06 (-5 ) 

1 . 20 (- 3 ) 

(8,1) 

1.61  (-5) 

1.43 (-5) 

1 . 35  (- 5 ) 

4.31  (-4) 

Diagonal 

Components 

(2,2) 

5 . 90  (-3) 

6.37  (-3) 

3 . 99 (-3) 

9. 55 (-3) 

(3,3) 

1 . 12  (-3) 

1.14  (-3) 

1.28  (-3) 

3.09 (-3) 

(4,4) 

1.321-4) 

3.69 (-4) 

4 . 05  (-4 ) 

2 . 35  (-  3) 

(6,6) 

5.38  (-6) 

7 . 00 (-5) 

9 . 30 (-5) 

1.37  (-3) 

(8,8) 

5.64  (-7) 

2.07 (-5) 

1 . 38 (-5) 

7.94  (-4) 

61 


Table  7.  Maximum  amplitude  in  z  of  the  Fourier  components  of 
the  x-velocity  u<m,n,z,t)  in  the  mixed  spectral 
(x-y)  -  physical  (z)  space  representation  for  Run  19. 


Time  (t) 


Component  (m,n) 

30 

45 

60 

75 

Two-Dimensional 

Harmonics 

(1,0) 

4.87  (-3) 

2.08  (-3 

1 . 40  (-3) 

2 . 58  (- 3) 

(2,0) 

8.51 (-4) 

2.57  (-4 ) 

1 . 4  2  (-3 ) 

1 . 85  (-3) 

(3,0) 

8.53  (-5) 

2.52  (-5) 

1 . 87 (-4 ) 

8 . 56  (-4 ) 

(4,0) 

1.92  (-5) 

2 . 4  2  (-6) 

7.07 (-5) 

2.98  (-4) 

(6,0) 

2 .70  (-6) 

3.29 (-7) 

1 . 20  (-6) 

4.39  (-5) 

Oblique 

Harmonics 

(1,1) 

2 . 76  (-2 ) 

6. 30 (-3) 

1.11 (-3) 

4.97  (-3) 

(2,1) 

1.31  (-3) 

9 . 32 (-5) 

3 . 95  (-4 ) 

2 . 96  (-3) 

(3,1) 

9.69  (-5) 

4.41 (-5) 

4.82  (-5) 

1 . 37  (-3) 

(4,1) 

2.29 (-5) 

3.91 (-6) 

3 . 30  (-5) 

3.92  (-4 ) 

(6,1) 

3.50 (-6) 

2.47 (-7) 

1.17 (-6) 

3 . 4  6  (-  5) 

Transverse 

Components 

(0,1) 

1 . 32  (-2) 

1 . 63  (-2 ) 

1.95  (-2) 

1.95  (-2) 

(0,2) 

2.57  (-2) 

4.19  (-2) 

4.08  (-2) 

3.65 (-2) 

(0,3) 

6.14  (-3) 

9.12  (-3) 

5.51  (-3) 

1.14 (-2) 

(0,4) 

1.91  (-3) 

6.34  (-3) 

5.62 (-3) 

7.72 (-3) 

(0,6) 

1.38  (-4 ) 

1.53  (-3) 

3 . 08  (-3 ) 

6 . 09 (- 3 ) 

Diagonal 

Components 

(2,2) 

1 . 30  (-3) 

3.55  (-4) 

2.44 (-3) 

3.55 (-3) 

(3,3) 

1.14  (-4) 

8.14  (-5) 

1 .45  (-4) 

1.62 (-3) 

(4,4) 

1 . 37 (-5 ) 

8.86  (-6) 

2.20  (-5) 

4 . 22 (-4 ) 
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Figure  Captions 


Fig.  1.  A  plot  of  the  profile  u(z)  of  the  x-velocity  of 
the  two-dimensional  initial  disturbance  used  in  Run  2.  This 
disturbance  is  chosen  as  the  least  stable  eigenmode  of  the 
Orr-Sommerfeld  equation  for  plane  Poiseuille  flow  with 
R  =  2935  and  a  =  1.3231.  The  phase  of  the  initial 
perturbation  is  chosen  so  that  the  maximum  velocity 
perturbation  occurs  initially  at  x  =  0,  where  this  plot 
is  made.  The  initial  disturbances  used  in  Runs  3  and  4  are 
proportional  to  that  imposed  in  Run  2. 

Fig.  2.  A  plot  of  the  time  evolution  of  the  maximum  amplitude 
in  z  of  the  x-velocity  component  of  the  two-dimensional 
primary  disturbance  [that  depending  on  x  like  exp(iax)] 
and  its  harmonic  [depending  on  x  like  exp(2iax)]  for 
Run  2.  At  late  times  the  primary  disturbance  is  undergoing 
slow  growth. 

Fig.  3.  A  plot  of  the  z-profile  of  the  x-velocity  component 
of  the  primary  disturbance  in  Run  2  at  t  =  120. 

Fig.  4.  A  plot  of  the  mean  velocity  u(z)  in  Run  2  at  t  =  120. 

For  comparison,  the  undisturbed  plane  Poiseuille  flow  profile 
2 

1-z  is  also  plotted. 

Fig.  5.  A  plot  of  the  curvature  -u,((z)  of  the  mean 

velocity  profile  in  Run  2  at  t  =  120.  For  comparison,  the 
curvature  2  of  the  parabolic  profile  1-z  is  also  plotted. 

Fig.  6.  A  plot  of  the  maximum  amplitude  A  of  the  primary 
disturbance  and  its  harmonic  for  Run  4.  See  Fig.  2. 

Fig.  7.  A  plot  of  the  maximum  amplitude  A  of  the  primary 
disturbance  and  its  harmonic  for  Run  1  at  R  =  2500.  The 
initial  disturbance  is  chosen  to  be  the  least  stable  eigen- 
mode  of  the  Orr-Sommerfeld  equation  with  R  =  2500  a  -  1.3231  . 

See  Fig.  2. 
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Fig.  8.  A  plot  of  the  z-profile  of  the  x-velocity  component 
of  the  two-dimensional  primary  disturbance  u(z)  for  the 
plane  Couette  flow  Run  6  at  R  =  5000,  a  =  2.  The  initial 
conditions  are  chosen  as  the  symmetrized  combination  (3.2) 
of  the  least  stable  eigenmode  of  the  Orr-Sommerf eld  equation 
and  its  complex- conjugate  reflected  eigenmode. 

Fig.  9.  A  plot  of  the  maximum  amplitude  A  vs  t  of  the 
x-velocity  of  the  primary  disturbance  and  its  harmonic  for 
the  two-dimensional  plane  Couette  flow  Run  6  at  R  =  5000. 

Fig.  10.  A  plot  of  the  maximum  amplitude  A  vs  t  for  the 
two-dimensional  primary  disturbance  and  its  harmonic  and 
the  three-dimensional  primary  wave  [witn  a  =  1.3231  and 
8=1]  in  Run  3A.  The  two-dimensional  part  of  the  initial 
conditions  is  the  same  as  in  Run  3  of  Sec.  3  while  the  three- 
dimensional  component  is  introduced  initially  through  round¬ 
off  error.  Note  that  the  amplitude  of  the  three-dimensional 

O 

component  is  multiplied  by  10  . 

Fig.  11.  A  plot  of  the  z-profile  of  the  x-velocity  of  the 
two-dimensional  disturbance  used  in  Runs  7-13.  This  disturb¬ 
ance  is  chosen  to  be  the  least  stable  eigenmode  of  the  Orr- 
Sommerfeld  equation  for  plane  Poiseuille  flow  at  R  =  1250 
with  a  =  1 .  The  amplitude  of  this  initial  perturbation  is 
reduced  by  25%  in  Run  9. 

Fig.  12.  A  plot  of  the  z-profile  of  the  x-velocity  of  the  initial 
three-dimensional  disturbance  imposed  in  Runs  7-9.  This  dist¬ 
urbance  is  chosen  as  the  least  stable  eigenmode  of  the 
Orr-Sommerfeld  equation  at  R  =  1250  with  a  =  1  and  8=1. 

In  Run  9,  the  initial  amplitude  is  decreased  by  25%. 

Fig.  13.  A  plot  of  the  maximum  disturbance  amplitudes  A  vs 
t  for  the  three-dimensional  plane  Poiseuille  flow  Run  7  at 
R  =  1250.  The  plotted  amplitudes  are  the  maxima  in  z  of 
x-velocity  of  the  two-dimensional  primary  disturbance  and  its 
harmonic  (multiplied  by  10)  and  the  three-dimensional  primary 
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disturbance  [depending  on  x  like  exp { iotx+iBy)  with 
a  =  1  and  3=1]. 

Fig.  14.  Same  as  Fig.  13  except  for  Run  8.  The  amplitude 
of  the  two-dimensional  harmonic  is  not  multiplied  by  the 
factor  10  used  in  Fig.  13. 

Fig.  15.  Plots  of  the  mean  velocity  profiles  (4.2)  for 

Run  8.  _  t  =  30; -  t  =  45;  -  t  =  60; 

.  t  =  75  .  The  corresponding  plots  for  t  =  30  and 

t  =  45  for  Run  7  are  indistinguishable  from  those  plotted 
in  this  Figure. 

Fig.  16.  Plots  of  the  instantaneous  x-velocity  profiles  u (x,y,z,t)  vs. 

z  at  t  =  30  for  Run  7 .  -  x  =  0 ,  y  =  -  ir  ; - 

x  =  0/  y  =  —  17 / 2  ;  — - —  —  x  =  0,  y  =  0;  •»»»•»»»  * 

x  =  1 3  V 8  r  y  =  -5V8. 

Fig.  17.  Same  as  Fig.  16,  except  at  t  =  45.  Labeled  curves 
are  plotted  at  the  same  values  of  x  and  y  as  in  Fig.  16. 

Fig.  18.  Contour  plots  of  the  instantaneous  x-velocity  u(x,y,z,t) 
in  the  y-z  plane  at  x  =  0  in  Run  7.  (a)  t  =  0,  contours 
0,0. 9(0.1).  (b)  t  =  15,  contours  0,l(.l).  (c)  t  =  30, 

contours  0, 1(0.1).  (d)  t  =  45,  contours  0, 1(0.1).  Note  that 

while  the  tick  marks  along  the  y-axis  do  indicate  the  available 
y-resolution,  the  tick  marks  along  the  z-axis  are  included  for 
reference  only  and  do  not  indicate  the  available  resolution  of 
the  Chebyshev  series,  especially  near  the  walls. 

Fig.  19.  Contour  plots  of  the  x-vorticity  component  (x,y,z,t) 
in  the  y-z  plane  at  x  =  0  in  Run  7.  (a)  t  =  0,  contours 

-0.35,  0.35(0.07).  (b)  t  =  15,  contours  -0.32,0.32(0.08). 

(c)  t  =  30,  contours  -1,1(0. 2).  (d)  t  =  45,  contours 

-2, 2(0. 4) . 

Fig.  20.  Contour  plots  of  flow  components  in  Run  7  at  t  =  30 
in  the  y-z  plane  at  x  =  0.  (a)  v(x,y,z),  contours  -0.06,0.06 

(0.01).  (b)  w  (x,y,z),  contours  -0.024,0.012(0.004).  (c) 
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10  2  '  y '  z  ^  '  c°htour  s  -2. 5, 2(0. 5).  (d)  ^  ,y  ,z)  ,  contours 

-0.3,0.3(0.06). 

Fig.  21.  Contour  plots  of  the  instantaneous  x-ve]ocity 
u(x,y,z,t)  in  the  x-y  plane  at  z  =  cos  3n/16  in  Runs  7,8. 
(a)  t  =  0,  contours  0.19,  0.41(0.01)  .  (b)  t  =  30, 

contours  0.21,  0.42(0.01)  for  Run  7.  (c)  t  =  60, 

contours  0.16,  0.84(0.04)  for  Run  8. 

Fig.  22.  A  plot  of  the  z-profile  of  the  two-dimensional 
primary  disturbance  u(l,0,z,t)  [see  (2.3)]  in  Run  7: 

(a)  t  =  30;  (b)  t  =  45. 

Fig.  23.  A  plot  of  the  x-profile  of  the  two-dimensional 
harmonic  disturbance  u  (-2- ,  0 ,-ir,  t )  [that  depends  on  x  and  y 
like  exp(2iax)]  at  t  =  45  in  Run  7. 

Fig.  24.  A  plot  of  the  z-profile  of  the  primary  three-dim¬ 
ensional  disturbance  u(l,l,z,t)  [that  depends  on  x  and 
y  like  exp (iax+i8y) J  at  t  =  45  in  Run  7. 

Fig.  25.  A  plot  of  the  maximum  disturbance  amplitudes  vs  t 
in  Run  9  at  R  =  125U.  This  run  is  the  same  as  Run  8  except 
that  the  initial  amplitudes  of  both  the  two-  and  three- 
dimensional  disturbances  are  reduced  by  25%. 

Fig.  26.  A  plot  of  the  z-profile  of  the  least  stable  eigen¬ 
mode  of  the  Orr-Sommerfeld  equation  for  plane  Poiseuille 
flow  at  R  =  1250  with  a  =  1  and  8=4.  Observe  that 

this  mode  is  concentrated  near  the  center  of  the  channel, 
so  it  has  little  opportunity  to  interact  with  the  least 
stable  two-dimensional  disturbance  which  is  concentrated 
near  the  walls.  This  kind  of  three-dimensional  disturbance 
typically  does  not  lead  to  transition.  The  eigenvalue  of 
the*Orr-Sommerfeld  equation  is  w  =  0.9076  -  i0.0579.  The 
phase  velocity  of  this  mode  is  about  0.9,  in  contrast  to 
the  phase  velocity  of  wall  modes  which  are  usually  in  the 
range  0. 3-0.4. 

Fig.  27.  A  plot  of  the  z-profile  of  the  second  least  stable 
eigenmode  of  the  Orr-Sommerfeld  equation  for  plane  Poiseuille 
flow  at  R  =  1250  with  a  =  1  and  6=4.  This  wall  mode 
has  eigenvalue  w  =  0.4635  -  i0.1607.  This  mode  is  used 
as  the  three-dimensional  disturbance  in  Run  13. 
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Fig.  28.  A  plot  of  the  maximum  disturbance  amplitudes  in  Run 
14  at  R  =  750.  This  run  uses  8  x  8  x33  spatial  resolution. 

Fig.  29.  Same  as  Fig.  28,  except  for  Run  15.  The  amplitude  of 
the  two-dimensional  harmonic  components  is  multiplied  by  a 
factor  10.  Run  15  differs  from  Run  14  only  in  that  the  spatial 
resolution  is  16  x  16  x  33.  Observe  that  the  flow  no  longer 
undergoes  'transition' . 

Fig.  30.  A  plot  of  the  maximum  disturbance  amplitudes  vs  t 
for  Run  16  at  R  =  500.  The  initial  conditions  for  this  run 
are  the  output  conditions  from  Run  8.  These  initial  conditions 
are  intended  to  simulate  turbulence. 

Fig.  31.  A  plot  of  the  z -profile  of  the  initial  two-dimensional 
disturbance  used  in  Run  17  and  18.  This  disturbance  is  const¬ 
ructed  from  the  least  stable  eigenmode  of  the  Orr-Sommerfeld 
equation  for  plane  Couette  flow  at  R  =  1250  with  a  =  1 
and  8=0  by  symmetrizing  according  to  (3.2). 

Fig.  32.  Same  as  Fig.  31,  except  that  the  profile  of  the 
initial  three-dimensional  disturbance  for  Run  17  and  18  is 
plotted. 

Fig.  33.  A  plot  of  the  maximum  disturbance  amplitudes  A 
vs  t  in  Run  17  for  plane  Couette  flow  at  R  =  1250. 

Fig.  34.  Same  as  Fig.  33,  except  for  Run  18,  also  at  R  =  1250. 

Fig.  35.  A  plot  of  the  mean  velocity  profile  u(z)  at  t  =  60 
in  Run  17. 

Fig.  36.  Same  as  Fig.  33,  except  for  the  plane  Couette  flow 
Run  19  at  R  =  1000.  Symmetrized  initial  conditions  are  used. 

Fig.  37.  Same  as  Fig.  33,  except  for  the  plane  Couette  flow 
Run  20  at  R  =  1000.  In  this  case,  unsymmetrized  initial 
conditions  are  used.  Observe  that  the  high  frequency  time 
oscillations  in  the  maximum  amplitude  of  the  two-dimensional 
primary  disturbance  and  its  harmonic  disappear. 
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